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1  Introduction 


Interest  in  nuclear  power  for  both  space  exploration  and  terrestrial  use  has  increased 
in  the  technical  community  in  the  recent  years.  Potential  space  missions  that  will 
demand  dependable  power  levels  beyond  conventional  sources  are  now  being  con¬ 
sidered  with  new  enthusiasm.  President  Bush  has  recently  expressed  his  desires  for 
the  nation  to  become  committed  to  putting  a  space  station  in  orbit  and  developing 
a  lunar  ba^  that  will  act  as  stage  for  the  exploration  of  Mars.  There  are  sugges¬ 
tions  within  the  civilian  community  that  the  development  of  manufacturing  plants 
in  space  is  receiving  strong  consideration.  Nuclear  power  is  seen  m  a  feasible  source 
of  the  energy  needed  to  accomplish  these  goals. 

Likewise,  developing  strong  concerns  of  possible  greenhouse  effects  from  using 
fossil  fuels  is  renewing  interest  in  developing  small  terrestrial  nuclear  reactors  with 
passive  safety  designs.  Even  some  established  critics  of  nuclear  power  suggest  that 
nuclear  power  would  be  acceptable  if  reactor  designs  with  new  passive  safety  features 
are  developed.  It  is  evident  that  nuclear  reactors  will  be  included  in  future  energy 
considerations. 

Therefore,  development  of  new  analytic  tools  is  desirable  to  enhance  reactor  de¬ 
signs.  This  research  resulted  in  developing  a  new  finite  element  code  which  will  be 
used  in  reactor  design,  development,  and  analysis,  or  for  any  time  dependent  radia¬ 
tion  transport  problem  defined  in  XZ  or  RZ  geometries. 

1.1  Research  Objectives 

The  object  cf  tliis  research  was  the  development  of  a  new  neutron  and  gamma  ra¬ 
diation  transport  code  with  the  followin^j  characteristics.  It  must  be  able  solve  two 


dimensional  time  dependent  problems  in  either  XZ  slab  or  RZ  cylindrical  geometry. 
The  angular  flux  must  be  defined  so  that  some  anisotropic  flux  behavior  can  be 
modeled.  The  spatial  discretization  must  be  adaptable  to  both  XZ  and  RZ  geome¬ 
tries,  as  well  as  future  options  that  may  be  incorporated  later.  The  code  must  have 
the  ability  to  model  systems  that  have  strong  energy  dependency.  And  finally,  if  a 
system  is  modeled  involving  one  or  more  fissile  isotopes,  the  source  must  include  a 
delayed  neutron  contribution. 

The  code  will  be  used  to  investigate  the  following: 

a.  Neutron  wave  effects. 

b.  Spatially  dependent  subcritical  and  critical  source  driven 
transients. 

c.  The  effect  on  accuracy  of  tradition  approximations  in 
solving  time  dependent  transport  problems.  This  will 
include  comparisons  with: 

i.  The  diffusion  approximation. 

ii.  Treatment  of  the  fission  source. 

iii.  Treatment  of  upscatter. 

The  code  will  then  be  benchmarked  and  some  few-group  problems  identified  and 
solved  that  demonstrate  the  code’s  validity  and  applicability. 

FMP2DT  (Finite  element.  Multigroup,  P„,  2-Dimensional,  Time  dependent)  was 
developed  to  met  these  objectives.  It  solves  time  dependent  neutron  and/or  gamma 
radiation  transport  problems  in  XZ  and  RZ  geometries.  Finite  elements  were  se¬ 
lected  to  implement  spatial  discretization,  and  the  time  discretization  was  done  using 
Euler’s  backward  differencing  scheme. 

The  finite  element  scheme  is  employed  by  expanding  the  solution  of  a  set  of  local 
partial  differential  equations  with  a  set  of  bcisis,  or  interpolating  functions.  Using  the 
Galerkin  procedure,  interpolating  and  weighting  functions  were  chosen  that  adapt 
to  both  of  the  stated  geometries,  or  any  geometry  option  added  later,  no  matter  how 
irregular  it  may  be.  The  finite  element  method  has  a  firm  theoretical  foundation 
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which  guarantees  convergence  of  the  approximate  solution  [Ref.  4]. 

Euler’s  backward  differencing  method  was  used  for  the  time  discretization.  This 
is  an  implicit  scheme  that  is  numerically  stable  for  any  time  step  [Ref.3]  [Ref.  26]. 

The  angular  dependency  of  the  angular  neutron  flux  was  modeled  using  spherical 
harmonics.  Spherical  harmonics  form  a  complete  set  of  functions  that  describe  the 
angular  dependency  of  the  neutron  direction  [Ref.  1].  Spherical  harmonics  yield 
solution  results  of  arbitrarily  high  degrees  of  accuracy  depending  on  the  expenditure 
of  labor  to  do  the  resulting  calculations  [Ref.  2).  At  least  a  Pi  approximation  is  a 
necessary  requirement  to  observe  neutron  wave  behavior  and  model  some  anisotropic 
behavior  of  the  flux.  All  of  this  applies  to  the  gamma  flux  also.  Therefore  it  was 
modeled  likewise. 

1.2  Literature  Search 

A  literature  search  conducted  in  July  1989  revealed  no  previous  finite  element  neu¬ 
tron  code  developments  with  delayed  neutron  sources  and  inhomogeneous  sources. 
There  are  numerous  codes  that  model  time  dependency  with  other  discretization 
characteristics. 

Kinetics  codes  are  available  with  numerous  finite  differencing  schemes.  Monte 
Carlo  codes  are  available  in  which  a  detailed  spatial  model  of  the  reactor  can  be 
accomplished.  .Ml  of  these  codes  have  legitimate  applications  where  they  have  ad¬ 
vantages  over  other  methods.  They  also  have  their  constraints.  Monte  Carlo  codes 
require  a  statistical  approach,  and  experience  is  needed  to  ensure  validity  for  differ¬ 
ent  reactor  configurations.  The  finite  differencing  scheme  is  the  most  used  solver. 
For  stiff  problems,  mesh  spacing  can  be  complicated,  and  implementing  various 
boundary  conditions  can  be  tedious. 

The  finite  element  code  developed  by  this  research  is  more  flexible.  The  ap- 
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proximating,  or  interpolating,  functions  used  in  the  code  allows  the  incorporation  of 
complicated  geometries.  Boundary  conditions  are  carried  with  each  finite  element, 
and  therefore  somewhat  easier  to  implement  than  for  other  schemes.  In  fact,  imple¬ 
mentation  of  the  boundary  conditions  appears  to  be  the  most  attractive  feature  of 
finite  element  solution  schemes. 

Finite  element  codes  have  been  used  in  neutron  transport  codes  for  sometime. 
Several  codes  exist  for  steady  state  analysis.  FEMPID  (Finite  Element  Multi-group 
P„  1-DimensionaI)  is  a  radiation  transport  code  for  infinite  slab  geometry.  Buckling 
height  corrections  are  needed  to  adjust  for  leakage.  This  code  is  very  cost  effective  . 
FEMP2D  is  a  two  dimensional  version  which  analyzes  a  steady  state  configuration. 
It  is  a  Pi  code  designed  to  handle  XZ,  RZ  and  R0  geometries.  Likewise,  FEMP3D 
is  three  dimensional.  These  codes  are  all  written  for  vector  machines  and  are  coded 
in  FORTRAN  77. 

PERT2D  is  a  finite  clement  perturbation  code  that  models  small  changes  in  re¬ 
activity  [Ref.  19].  A  one  dimensional,  time  dependent,  finite  clement  a)dc,  TDFID, 
is  also  available,  but  it  does  not  have  a  delaj'cd  neutron  precursor  source  contri¬ 
bution  [Ref.  25].  SIILDTEMP  is  a  coupled  radiation  transport  and  temperature 
distribution  code  written  at  The  University  of  New  Mexico  [Ref.  24). 

J.  K.  Flctdier  has  suggested  finite  clement  options  in  some  of  his  transport  codes 
for  steady  state  [Ref.  11).  Finite  elements  have  been  used  to  discretize  the  angular 
dependency  of  the  neutron  direction  in  other  neutron  codes.  So  the  foundation  for 
using  finite  elements  in  neutron  transport  is  strong.  However,  the  literature  seardr 
did  not  reveal  any  previous  finite  element,  spatially  discretized,  multidimensional 
transient  codes.  Therefore,  the  development  of  this  new  code  had  a  strong  theoret¬ 
ical  foundation,  and  it  contributes  to  the  work  previously  done. 
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1.3  Importance 


The  importan.e  of  this  effort  is  that  a  finite  element  code  now  exist  for  time  depen¬ 
dent  radiation  transport  and  nuclear  reactor  analysis.  The  code  is  both  relatively 
easy  to  set  up  and  economical  to  use.  Benchmarking  shows  that  it  achieves  an  ac¬ 
ceptable  degree  of  accuracy  for  XZ  and  RZ  geometries.  Delayed  neutrons  can  now 
be  modeled  in  the  source  terms.  Although  delayed  neutrons  have  an  extremely  small 
population  in  reactors  compared  to  prompt  fission  neutrons,  their  presence  ensures 
that  the  reactor  is  controllable.  Thus,  developing  a  finite  element  code  that  con¬ 
siders  their  source  contribution  is  a  credible  enhancement  and  contribution  to  the 
inventory  of  codes  now  available.  This  code  will  make  a  contribution  to  fundamental 
engineering  knowledge.  It  is  also  a  valuable  stepping  stone  for  the  development  of 
higher  order  approximations  in  future  research. 
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2  The  Time  Dependent  Transport  Equation 


Neutron  transport  theory  y’-  .s  a  mathematical  expression  which  describes  neutron 
interactions  in  a  given  medium.  For  any  arbitrary  volume,  the  time  dependent 
neutron  transport  equation  may  be  written  in  the  following  form: 

-  ^  $(r,  E,  ^)  +  n .  V  $(r,  E,  Q,  t)  +  St(r,  E)  $(r,  E,  Q,  t) 

=  j  j  T.,{v,E'  E,^'  -^n)^T,E',€l',t)d9! dE' 

E' 

+  ^^(l-;5)  J  J  E')^r,E', a', t) da' dE' 

E' 

+  S  E  xm  X,  ft(r,  t)  +  S,„(r,  E,  S2,  i)  (1) 

Where  Equation  (1)  states  the  relationships  between  source  and  loss  terms  such 
that: 


{The  time  rate  of  change  of 
the  neutron  angular  flux 

neutron  1  i  /  prompt  fission  1  i  f  delayed  1 
inscattering  /  \  neutrons  J  1  neutrons  j 

,  f  inhomogeneous  1 
\  sources  J  ‘ 

The  streaming  and  total  interaction  variables  represent  loss  terms,  where  the  terms 
on  the  right  hand  side  represent  sources.  The  angular  neutron  flux,  $(r,E,0,<), 
is  a  function  of  spatial  variables  incorporated  in  r,  neutron  energy  E,  angular  vari¬ 
ables  incorporated  in  12,  all  of  which  are  evaluated  at  time  t.  12  is  represented  by 
an  azimuthal  angle  (j)  and  a  polar  angle  0  in  an  orthogonal  coordinate  system.  The 
polar  angle  is  usually  described  in  terms  of  its  cosine.  A  complete  set  of  functions 
that  describe  the  angular  dependency  of  the  neutron  flux  are  spherical  harmon¬ 
ics  [Ref.  1].  For  XZ  and  RZ  geometry,  with  azimuthal  symmetry,  the  angular  flux 


},  f  neutron  1  i  f  total  1 

\  streaming  J  \  interactions  J 


may  be  completely  described  with  the  spherical  harmonic  expansion  as, 

°o  ^  01  A-  1 

— — Pi^(cos6)il)im{r,E,t)  cos(m^)  (2) 

t=0  m=Q 

Likewise,  the  inhomogeneous  source  may  be  expanded  as, 

°°  '  9/  -I- 1 

5;„(r,  E,n,t)  =  Y^  ^  — —  Pr{cos  e)  Simir,  E,  t)  cos(m^)  (3) 

1=0  m=0 

Pl^{cosO)  are  associated  Legendre  functions  of  degree  /  and  order  m  [Ref.  5].  They 
are  orthogoi:?!  on  the  interval  ^  =  0,  to  0  =  it  radians.  They  obey  the  recurrence 
relationships  [Ref.  11]: 

(21  +  1)  cos OPricosO)  =  (/-m  +  l)P,:;!i(cos0)  +  (/  +  m)P,’:!i(cos<>)  (4) 
{2l  +  l)smOPr{cosO)  =  P{lt\cos0)-P;2t\cosO)  (5) 

(2/  +  1)  sin  0 P,'”(cos 0)  =  (/  +  m)  {/  +  m  —  1)  P,^7^(cos  0) 

-  (/  -  m  +  1)  (/  -  m  +  2)  P,!?r'(cos  0)  (6) 


The  streaming  term  in  Equation  (1)  is  •  V$(r,P,  This  term  is  a  rep¬ 
resentation  of  the  rate  of  change  of  the  neutron’s  angular  flux  along  a  streaming 
path  [Ref.  12].  The  •  V  term  is  sensitive  to  the  geometry  in  question.  For  XZ 
geometry,  it  is  defined  as, 

.  V  =  cos<?— -b  sin0  cos^—  (7) 

For  RZ  geometry,  with  symmetry  in  its  azimuthal  coordinate,  it  is  defined  as, 

r-,  «  5  ,  5  sin  sin  5 

.  V  =  cos  —  -b  sin  cos  ^  :;r - (8) 

az  or  r  dtp 

The  inscattering  probability  (or  cross  section)  in  Equation  (1)  may  be  represented 


9/4-1 

->  P,f2'  fi)  =  X;  -7^  S„(r,P;'  E)  P,(f2'  •  (9) 
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By  the  addition  theorem  [Ref.  6],  we  may  write 


PliP,'  ■  fl)  =  Pi(cos6)Pi{cos0')  +  2  V  II  P,”*(cos5)P/"(cosO')  cos[m(.^ -  ^')1  (10) 

[I  +  Jh)! 

m=:l 

Therefore  the  ihscatter  cross  section  may  be  expressed  as, 
s,(r,  = 

|p,(cos(?)P,(cos0')  +  2  Y,  7J^^r(“s®)^’r(cos5')  cos[m(^ - ^')]|  (H) 

(=0  V  m=l  J 

Exact  solutions  to  Equation  (1)  are  possible  in  only  a  few  special  cases.  In  prac¬ 
tice,  approximations  are  made  to  the  transport  equation  to  generate  solutions  that 
are  accurate  enough  for  specific  physical  interaction  processes.  We  now  consider 
the  development  of  an  approximate  solution  to  the  transport  equation  for  XZ  slab 
geometry. 
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3  Solution  in  XZ  Geometry 


To  get  an  approximate  solution  to  the  transport  equation  in  XZ  geometry,  the 
correct  substitutions  shown  above  are  implemented  into  Equation  (1).  The  transport 
equation  is  then  transformed  into  the  following  form: 

X  X  {cosO)i)imcos{m<f>) 

^  1=0  m=0 

f)  °°  '  2/4-1 

+  COS  0  ^  — —  P,”*(cos  6)  cos(m^) 

d  °°  *  2/  -b  1 

-b  sin  0  cos  ^ ^  Y  “7 — P/"(cos/?)  ■0/^  cos(7n^) 

/=0  m=0 

CO  I  2/4-1 

+  £  Zy  “7 — Pr(cos//)^/,„cos(m^) 

/=0  m=0 

r  r  97  4- 1 

=  //  i:-^S„P,(cos/?)P,(cos/?') 

JS'  J2'  '=° 


+ 


CO  I  2/  -I-  1 

■Y  E  ^^Pri^osO')^imCos{mcl>')da'dE' 

1=0  m=0 

/  /  E  ^  (2)  £  f,"(cos  «)  P,”(«>s  O')  cos[m{#  -  ^')) 


S'  ' 

oo  / 


•EE 

1=0  m=0 


21  +  1 
47r 


P;'"(cos  O')  'tjfim.  cos{m<f>')  do!  dE' 


+  ^  (1  — /?)  I  f  EE^  P;’"(cos  //')  cos{ni<f>')  dfi'  dE' 

^  ^  1=0  m=0 


S'  J2' 

~  '  2H-1 


+  4. 

/=0  m=0 


P;”*(C0S//)  S'/,„COS^  - 


+  (12) 

k=l 

Equation  (12)  can  be  imn^ediately  simplified  by  defining  fission  as  an  isotropic 
event.  This  has  been  shown  to  be  experimentally  correct  [Ref.  7],  and  it  sets  /  = 
m  =  0  in  the  fission  term.  Some  of  the  €1'  terms  can  also  be  simplified  by  integrating 


over  all  directions.  This  integration  is  defined  as, 


0  2ir 


j  dQ' =  j  J  dtj/ d(cosO'):= 


(13) 

J  J 
TT  0 

and,  it  can  be  immediately  used  in  both  the  fission  and  inscattering  terms.  Asso¬ 
ciated  Legendre  functions  are  orthogonal.  Integration  of  these  functions  is  defined 

J  PT (cos 0) pH {::osO)d(cose]  =  (M) 

where  6ik  and  S^n  are  Kronecker  delta  functions.  Since  Pq{cosO)  =  1.0,  multiply¬ 
ing  any  term  in  Equation  (12)  by  Pq  (cos  0)  will  not  change  the  equation’s  value.  To 
implement  orthogonality  (at  a  later  time),  the  precursor  term  and  the  fission  source 
term  can  be  multiplied  by  Pq{cos9).  Using  the  above  information  yields, 

1  Q  OO  i 

“  ^  S  S  (2^  +  1)  -P/'”(cos  e)  ipim  cos{m(j>) 

^  1=0  m=0 


OO  / 


di>h 


+  ^  (2^  cosO  Pi^{cos  0)  cos(m^) 

f=0  m=0 

+  13  X3  (2^  +  1)  smO Pi^{cos0)  cos(m^)  cos  (f> 

1=0  m=0 

OO  / 

+  ^‘13  13  (2^  +  1)  Pr(cosfl)  jj)!^  cos(m^) 

/=0  m=0 

OO  OO 


*  w 

=  /  E(2'  +  l)2.,fl(cos(!)£*„iB' 

jjji  1=0  1=0 

+  /  Y,<?l  +  l)^.,Pr(cosO)  E  E  cos(mf)dB' 

j^i  l=l  /=!  ni=l 

+  X,{i-0)  J  isSf^odE'PSi  cos  0) 

E' 

OO  I 

+  E  E{2'  +  l)^',”(“S»)S,„cos{m#) 

/=0  m=0 

+  Y,XkhGkP^[cosO) 

k=\ 


(15) 


Note:  Pl^{cos0)  =  Pi(cos9)  if  m  =  0.  The  following  trigonometric  identity  may 
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now  be  implemented: 


cos(m^)  cos(f>  =  ^  cos(m  +  1)^  +  ^  cos(m  —  1)^ 


(16) 


We  can  substitute  Equation  (16)  and  implement  recurrence  relations  defined  by 
Equations  (4),  (5),  and  (6)  and  manipulate  so  Equation  (15)  becomes, 

0  A  CO  / 

“  £  £  (2^  +  1)  -Pr(cos  9)  i^im  cos(m^) 

+  (2)  ]££)(/-  m  +  1)  P,^i(cos  0)  cos(m9i) 

+  (2)  £  £  {I  +  m)  PiU-^icqs  9)  cos(m^) 

;=0  m=0 

+  £  £  ■f/+i'^(cos  <?)  cos(m  +  1)^ 

i=0  rn=0 

"  £  £  Pi!lt\cos9)  cos(m  +  1)^ 

(=0  m=0 

+  £  £i  (^  +  m)(/  +  m  -  1)  P/!!!7^(cos  0)  cos(m  -  1)^ 

/=0  m=0 

-  X)  ]£  (^  “  "^  +  -  m  +  2)  Pm  f ^(cos  0)  cos(m  -  1)^ 

/=0  m=0 

00  / 

+  (2)  S,  E  E  (2'  +  1)  J'rCcos  0)  i>,„  cos(m^) 

/=0  m=0 

/OO  00 

E(2;  +  l)S„fl(cosO)E*o<iB' 

/=0  /=0 

/OO  00  / 

X:(2/  +  l)  S„Pr(cos<?)  £  cos{m^)dE' 

£!  1=1  1=1  m=l 

+  {2)Xp(l-;3)  J  i>^,<PoadE' PiicosS) 

E' 

00  I 

+  (2)  £  £  (2/  + 1)  Pr(cos  ^?)  Sim  cos(m^) 

/=0  m=0 

+  (2)i:XA:AfcC'fcPoVs<?)  (17) 

k=l 

Next,  the  subscripts  and  superscripts  may  be  set  to  put  the  transport  equation  back 
in  terms  of  spherical  harmonics.  This  can  be  done  because  starting  an  infinite  series 
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expression  at  a  different  point  will  not  change  series  c  >nvergence  as  long  as  all  the 
series  indices  are  likewise  changed.  Accomplishing  this  ])roduced  this  form: 

"  ^  £  S  (2^  +  1)  i^r(cos  0)  cos(m^) 

+  (2)  £  Z)  (/  -  rn)  (cos  9)  cos(m^) 

/=1  m=0 

+  (2)  Z  Z  -  cos(m^) 

;=0  r7i=0 

+  Z  Z  ■Pr(cos{/)  2^^^^^  cos(m(^) 

-  Z  Z  J^r(cosil)  cos(m^) 

Z=-l  m=l 

+  Z  Z  +  2)(^  +  m  +  1^ P,”*(cos cos(m^) 

Z=0  m=0 

—  ^  ^  (/  —  m  —  !)(/  —  m)  PI^(cos  0)  cos(m^) 

.’=1  m=-l 

oo  Z 

■h  (2)  St  2  (2^  +  1)  -Pz"‘(cos  <?)  V’Zm  cos(m<^) 

/=0  m=0 

/oo  CO 

J^{2l  +  l)S„Pi{cosO)Y,i>lodB’ 

i~o  /=o 

/oo  oo  I 

i:(2(  +  l)  S„P,’”(cos«)  i:  £  cos{m4)dE' 

j^i  Z=1  Z=1  ni=l 

+  {2)Xp(l-«  J  uEj^aodE' {^(msO) 

E' 

00  Z 

+  (2)  E  E  (2(  +  l)P,’”(cos«)S,„co8(m«l) 

Z=0  m=0 

+  (2)ExtAiftP»(cosO)  (18) 

Zi:=l 

Equation  (18)  is  still  in  exact  form.  However,  the  summations  to  infinity  prohibit 
practical  implementation.  The  expansion  coefficients  0zm  are  what  we  are  ultimately 
solving  for.  They  can  be  found  by  integrating  over  all  directions  and  using  orthog¬ 
onality.  If  every  term  in  Equation  (18)  is  multiplied  by  P^{cosO)  cos{N(f>),  and  the 
integration  over  all  directions  is  done,  there  are  two  orthogonal  relationships  to  ad- 
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dress.  One  has  been  defined  by  Equation  (14).  The  other  is  related  to  the  azimuthal 
angle,  (j).  It  is  defined  as: 

r  0  if  m  ^  N 


2‘7r 

J  cos(m^)  cos{N <f>)  d(f>  = 


TT  if  m=-  N 

27r  if  m  =  iV  =  0 


(19) 


Now  we  multiply  Equation  (18)  by  Pj^  {cos  6)  cos{N(f>),  and  set  it  up  to  integrate 
over  all  directions.  This  yields: 

2  d  ^  ' 

/  /  E(2/+i)  P;’”(cOS  6)  cos(m^)  P^{cos  $)  cos{N4>)  d^  d{cos  0) 

J  V  at 


0  27r 


00  / 


+  (2)  /  E  E  -Pr(cos0)  cos{m(l>) 


1=1  m=0 

•P^(cos  0)  cos{N<f)  d(f)  d{cos  0) 


dz 


0-  27r 


oo  I 


+  (2)  /  /  E  E  +  1)  -PrCcos  /?)  — cos(m^) 


;  5  /=0m=0 

•Pjl (cos  0)  cos{N4>)  d<l)  d{cos  0) 


r  r  ‘ 

+  J  J  Y2  P/”(cos0) — cos(m^)P^(cos<?)  cos{N <f>)  d(f>  d{cos  0) 

It  0 

»  00  /  Q  ; 

“  /  /  E  E  P/"(cos<?) — cos{m(j))P^{cosO)  cos{N4')d^d{cos6) 
i  I  i=-i  m=i 

+  J  y^EE(^  +  ”^  +  2)(/  +  m  +  1)  P;’"(cos  fl)  cos{m^) 


•Pk{cosO)  cos{N4>)  d<j>d{cosO) 

-  j  j  H  Y.  (/-m-l)(/-m)P,'"(cos<?) 


TT  0 


/=!  m=-l 

P;^  (cos  <?)  cos{N<f>)  d<f>  d{cos  0) 

OO  / 


dfpl—lm+l 

dx 


cos(m^) 


0  2jr 


+  (2)  /  /  St  E  E  (2^  +  1)  Pr(cos  ^1)  ipim  cos{m<f>)Pj^ (cos  <?)  cos{N<j>)  d(f)  d{cos  0) 

i  Q  1=0  m=0 


0  2ir 


=  (2)  /  /  /  E  (2^  +  1)  S„  P/(cos  <1)  E  i’lo  dE'Pj^ (cos  0)  cos(IV^)  d(cos  tf) 

i  0  E'  '=0 


/=0 


0  2;r 


+  (2)  /  /  /  E(2'  +  l)S.,/’"(cosfl)£  2 

Tt  0  E'  m=l 
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‘P^{cosO)  cos{N (f>)  d<f>  d{cos  0) 

0  2ir 

+  (2)/  j  xAi-?) 

TT  0 

°  oo  I 

+  (2)  /  /  (2/  +  1)  PP{cos  0)  SimCOsim(j))P^ (cos  0)  cos{N(f>)  d(f>  d{cos  0) 

i  i  1=0  m=0 

0  2x  n 

+  (2)  j  j  ilXk  h  Ck  P^icos  0)P^ (cos  0)  cos(N(l>)  d<f>  d{cos  0)  (20) 

v  0 

So  far  no  approximations  to  the  transport  equation  have  been  made.  Simplifications, 
such  as  the  isotropic  nature  of  fission,  were  theoretically  correct.  However,  we  are 
now  ready  to  incorporate  some  approximations  to  the  transport  equation  to  get  it 
into  a  solvable  form. 

3.1  The  Pi  Approximation 

Equation  (20)  is  exact.  Now  the  first-  approximation  to  the  transport  equation  is  to 
be  made.  The  summation  over  I  is  truncated  to  some  finite  value  for  the  angular 
flux  and  inhomogeneous  source  expansions.  Truncating  the  upper  limit  of  /  also 
sets  the  upper  limit  of  m.  If  /  =  0,  then  m  =  0  in  all  cases.  If  /  =  1,  then  m  can 
take  on  values  of  0  and  1.  For  certain  values  of  /  and  m  some  summation  terms 
in  Equation  (20)  will  not  ex‘d«.  Making  the  Pi  approximation  states  that  either 
all  expansion  coefficients  0  or  that  the  partials  of  ^2m  in  either  the  x  ot  z 

directions  are  approximately  zero.  Therefore  this  approximation  ’s  not  uniquely 
defined. 

The  Pi  approximation  states  that  the  angular  flux  can  be  adequately  defined  by 

^(r,i?,n,t)  =  ^  Po  (cos  f?)^oo(r,  .£,<)  + :;^Pi°(cos  0)^1  o(r,P,t) 

4r7r  47r 

3 

+  - —  P/(cos  0)  i{ji  1  (r,  E,  t)  cos  ^ 

QilT 


j  uTijilfoodE'  Pq  {cos  0)PI(  {cos  0)  cos{N (j>)  d(f>  d{cos  0) 
E' 
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Some  values  for  the  associated  Legendre  polynomials  are  [Ref.  27]: 


Pq  (cos  9)  =  1  P°(cos  0)  =  cos  0 

Pl{cosO)  =  s\nO  (21) 


For  the  angular  flux  expansion  we  then  have, 
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$(r,  E,  t)  =  —  ^00  +  —  cos  0  ^1 0  +  —  sin  0  1  cos  ^  (22) 

47r  47r  47r 

where  the  function  arguments  of  the  expansion  coefficients  have  been  dropped.  The 
Pi  approximation  is  the  first  correction  to  the  diffusion  equation. 

Diffusion  theory  is  valid  in  large  homogeneous  or  nearly  homogeneous  reactors 
in  which  the  curvature  of  the  reactor  is  close  to  the  mean  free  paths  of  the  neu¬ 
trons  [Ref.  9].  Diffusion  theory  breaks  down  near  reactor  boundaries  or  strong 
absorbing  materials  [Ref.  10].  Therefore,  since  the  P\  approximation  resembles  dif¬ 
fusion  theory,  it  is  expected  that  similar  properties  would  hold  for  it.  However,  the 
diffusion  coefficient  should  be  better  defined  for  the  P\  transport  equation  approxi¬ 
mation. 

Now  if  we  let  if  =  AT  =  0  in  Equation  (20)  and  integrate  over  all  directions,  we 
get, 

~  ^  V’oo  +  — 1-  "b  ’/’oo  “  J  1^00  dE'  +  Xp{^  ~  P)  J  ^00  dE' 

E'  E' 

+  >Sbo  +  X/  (23) 

k=\ 

Doing  the  same  procedure  but  setting  K  =  \  and  iV  =  0  results  in. 


3  5  ,  d^o  ,  .  ov.  / 

U  ai  *  “ + ^ -ar + ^ -ar + ^  " 


3  J  S.j  0  dlS'  4"  3  5*1 0 

E' 


(24) 
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Applying  the  Pi  approximation  to  Equation  (24)  yields: 

^|v>io  +  ^  +  3S,^,„  =3  J  S„*„<(E'  +  3S,o  (25) 

E' 

Again,  setting  A”  =  A*  =  1  in  Equation  (20),  integrating,  and  implementing  the  Pi 
approximation  results  in, 

2|fe+^  +  3S,*,=3  /  S..A,rf£:'  +  3S.r  (26) 

B' 

(In  Equation  (26),  for  m  =  1,  the  ^oo  coefficient  is  multiplied  by  2.  This  is  caused  by 
not  using  the  recurrence  relationship  defined  by  Equation  (5)  since  for  the  case  when 
m  =  0,  Pi~^  {cos  0)  would  result  [Ref.  11].)  Three  equations  have  been  developed 
with  three  unknowns,  ^oo)  ^lO  and  ^n. 


3.2  The  Multigroup  Approximation 


This  approximation  entails  dividing  the  infinite  energy  spectrum  into  discrete  energy 
groups  that  are  defined  so  that  the  spatial  shape  of  the  flux  does  not  change  in 
the  discrete  groups  [Ref.  10].  If  the  energy  groups  are  defined  small  enough,  then 
this  isn’t  much  of  an  approximation  at  all.  The  multigroup  approximations  to 
Equations  (23),  (25),  and  (26)  yield: 


G  Tl 

3'=I  t=I 

f  ^A  +  ^  +  3Sf?(fo  =  3  E  + 

3  Q 


®  +  ^  +  3  s?  «.  =  3  E  Ej;-»  +  3  5f . 

3'=I 


Vg  dt 


(27) 

(28) 
(29) 


Each  of  the  above  equations  is  valid  for  a  particular  energy  group  g,  with  g  = 
1, 2,  •  •  • ,  G  as  possible  values. 
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3.3  Time  Discretization:  Euler’s  Backward  Differencing 


Equations  (27),  (28),  and  (29)  will  now  be  discretize  in  time  using  Euler’s  back¬ 
ward  differencing  scheme.  This  is  an  implicit  scheme  since  it  involves  variables  of 
the  present  time  step  on  both  sides  of  the  equation.  This  scheme  is  numerically 
stable  [Ref. 3],  which  is  its  most  desirable  characteristic.  The  present  time  step  is 
designated  with  a  (iV  +  1)  superscript,  and  the  previous  time  step  is  designated  with 
a  (N)  superscript.  Implementing, 

(W+l)  (N)  (W+l)  (iV+l) 

1  ql,  3  —11)3  fill,  3  fill)  ^  ^ 

Voo  Woo  ,  ovio  ovii  A.y3,j,  9  _  y'  y3'-9.u  s' 


Vn 


dz 


dx 


G  (w+i)  (w+i) 

+  4  (1  -  ^)  E  S?'  +  s„„’  +  E  X?  At  c{™'  (30) 

g'=\  k=l 


(iv+l)  (W)  (W+l) 

3  ^10  -Ao  ,  ® 


At 


‘'00 


dz 


+  32f>)’,o' 


(K+l) 

9 


(W+1) 


=  3  E  +35) 

9'=\ 


Itf+l) 

9 

10 


(31) 


(W+l)  (N) 

Vg  At 


+ 


(w+i) 
■(  3 
00 


dx 


-hSSf^ 


(A'+I) 
,  9 
11 


=  3E 


2,,. 


(W+1) 


(W+J) 

+  35,/ 


5'=1 


Now  all  values  with  a  superscript  (iV)  are  known  because  they  are  values  from  a 
previous  time  step.  The  ^  values  that  have  the  superscript  (A'+l)  are  our  unknowns. 
The  inhomogeneous  sources  are  always  known  no  matter  their  superscript  because  we 
assume  that  the  extraneous  neutron  source  centributions  can  be  directly  c..tculated. 

The  precursor  values  in  Equation  (30)  now  are  addressed.  The  precursor 

concentration  for  the  present  time  step  may  be  calculated  in  terms  of  the  previous 
time  step  precursor  concentration.  The  rate  of  change  of  the  precursor  concentration 
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with  respect  to  time  is  [Ref.  7], 


=  A  E *:  =  l,3.--,n  (33) 

a'=i 

Discretization  of  Equation  (33)  in  time  using  Euler’s  backward  differencing  yields, 

(34) 


Mi^’)  G  (A’+i) 

s'=i 


At 


Solving  Equation  (34)  for  yields, 

1  \  ^(N\  .  (  At 


= (tot)  + (tot) 


(.v+») 

s' 


(35) 


Now  the  present  time  step  value  of  the  precursor  concentration  can  always  be  calcu¬ 
lated  in  a  straightforward  manner  using  its  previous  time  step  value.  Substituting 
Equation  (35)  into  Equation  (30)  and  rearranging  yields, 


(A'+>)  (.v+i) 


+OT+  ax 


+  E?lfoo‘’ 


Uj  At 


+E 


1  CW>  G  (w+i) 

— Vo  +  E  Sr’l!-oo“  +Soo’ 

il’=l 

n  /  1  \  1  <"+'> 

.\|(i-/5)+EAiai  , ,  _L.  E  ■^E5'v''oo’' 

k=l  3^=1 

xfAtcr 


fcj  l+AiAi 


(36^ 


A  Taylor  series  expansion  can  show  that  the  error  for  this  backward  difference  ap¬ 
proximation  is  C?(Af)  [Ref.  37].  Thus  a  large  At  with  a  quickly  changing  process 
a)uld  yield  a  large  error.  However,  computational  cost  stipulate  the  selection  of  the 
largest  At  possible  that  yields  a  convergent  solution. 

If  there  arc  G  energy  groups,  then  G  equations  may  be  written  using  Equa¬ 
tion  (36).  Likewise,  the  same  may  be  done  using  Equations  (31)  and  (.32)  respec¬ 
tively.  Thus,  Equations  (36),  (31),  and  (32)  may  be  written  in  va:tor  form  repre- 
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senting  G  energy  groups.  The  vector  form  for  Equation  (36)  is, 


=  +  V  +  C(^)  (37^ 

here  and  C(^)  are  vectors  containing  G 

entries.  S°°,  V,  and  S-^  are  GxG  matrices.  All  these  terms,  and  other  vectors  and 
matrices  developed  later,  are  defined  in  Appendix  A.  Equation  (31)  in  vector  form 
is, 

^  +  3  S^®  =  3  y  'S'S?  +  3  (38) 

oz 

Equation  (32)  in  vector  form  is, 

4-  +  3  ^  =  3  V  +  3  (39) 

ox 

Now  solve  Equation  (38)  for 

3  S^®  =  3  y  +  3  -  ~ 

OZ 

Define  D^®  =  {S^®}  ^ ,  then, 

^  D'®  [3  y  +  3 -  4-  (40) 

o  [  Oz  \ 

For  Equation  (39),  solve  for  'Sf 

3  =  3  y  +  3  -  4- 

ox 

Define  D“  =  {S"}“\ 

=  \d^^  3 y 3 -  4- (41) 

0  Ox 

Define  =  4D^°  and  Then  substituting  Equations  (40)  and  (41) 

into  (37)  will  yield. 
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_  _^^(N+1)  .  y.00q^(iV+l) 

=  s^  +  s(^+')  +  y  9^J1'^  +  c(^) 


_ ^  nioy'Qpr(^) _ ^  niog(^+i) 

-|-i5“y*ff-|-Z)“S<r'* 

£/x  as 


(42) 


This  is  the  equation  in  which  the  spatial  discretization  using  the  finite  element 
method  will  be  made.  The  only  unknowns  are  in  the  vectors. 


3.4  Finite  Element  Implementation 

If  Equation  (42)  is  put  into  a  matrix  equation  form  of  say  Ax  =  b,  then  numeri¬ 
cal  algorithms  can  be  used  to  obtain  a  solution  [Ref.  18],  Thus  far,  the  transport 
equation  has  been  discretized  in  its  angular  dependency  using  spherical  harmonics, 
discretized  in  its  energy  dependency  using  the  multigroup  approximation,  and  dis¬ 
cretized  in  time  using  Euler’s  backward  differencing.  Now  we  turn  to  the  spatial 
discretization. 

The  finite  element  procedure  consists  of  approximating  a  solution  with  a  trial 
function.  The  set  of  functions  that  approximate  the  solution  vector  is  referred  to 
cis  a  trial  space.  Once  the  proper  trial  function  is  selected,  then  we  invoke  the 
method  of  weighted  residuals  incorporated  with  the  Galerkin  method.  This  will 
tend  to  spread  the  error  that  resulted  from  the  trial  function  approximation  so  that 
it  is,  in  some  sense,  small  over  the  whole  problem  domain.  For  transport  equation 
problems  such  as  this,  error  analysis  comparisons  between  the  finite  element  process 
and  finite  differencing  schemes  favor  using  finite  elements  [Ref.  19]  [Ref.  24].  We 
now  summarize  the  method. 


20 


Consider  the  following  equation, 

V2^-/  =  0  (43) 

where  /  is  a  known  function  of  the  independent  variables.  We  approximate  the 
solution  (j)  with  some  trial  function,  Most  likely  then  Equation  (43)  will  not  be 
true.  There  should  be  some  error  called  the  residual,  R,  such  that, 

VH-f  =  R.  (44) 

where 

M 

4  =  E<l>i  (45) 

t=i 

The  (j)iS  in  Equation  (45)  are  the  expansion  coefficients,  and  the  NiS  are  interpo¬ 
lating  functions  [Ref.  13].  If  the  residual  is  weighted  over  the  entire  domain  to  de¬ 
termine  the  s  such  that  the  error  is  small,  then  this  is  implementing  the  method 
of  weighted  residuals  [Ref.  13]  [Ref.  14].  We  now  choose  M  linearly  independent 
weighting  functions,  Wj,  so  that 

/  [V^ ^ - /]  WidD  =  J  RWidD  =  0,  (46) 

D  D 

where  D  in  Equation  (46)  denotes  the  problem  domain.  So,  in  some  sense,  the 
residual  i?  fs  0  over  the  entire  problem  domain.  The  Galerkin  method  simply  states 
that  the  weighting  functions,  W,,  may  been  the  same  as  the  approximating  or  trial 
functions  [Ref.  13]. 

Several  methods  exist  for  developing  the  finite  element  equations.  Using  the 
Galerkin  method  allows  the  development  of  these  equations  without  any  knowledge 
of  the  physical  processes  or  variational  calculus.  Imposing  it  allows  the  development 
of  a  numerical  algorithm  to  solve 

WidD  =  0.  (47) 

D 

The  Galerkin  method  is  now  used  to  solve  Equation  (42)  for  XZ  slab  geometry. 
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3.5  Finite  Element  Discretization  for  XZ  Geometry 


Linear  B  splines  will  be  used  for  the  interpolating  functions.  These  functions  are 
continuous  in  their  first  derivatives.  The  angular  flux  vectors  may  be  expanded  in 
theform, 

=  E  E  «!»”*  B,{x)  B,(z)  (48) 

p=l  9=1 

Ai  d  the  inhomogeneous  and  precursor  sources  may  likewise  be  expanded, 

b,(x)b,(z) 


p=l  9=1 

CW„cW  =  EEcOTB,(i)B,(z) 

P=1  9=1 


(49) 

(50) 


A  and  B  in  the  above  summations  represent  the  upper  limits  of  the  mesh  spacing 
in  the  x  and  z  directions  respectively.  Likewise,  the  p  and  q  subscripts  represent 
the  mesh  point  in  the  x  direction  and  the  q^^  m.esh  point  in  the  2  direction. 
Bp{x)  and  Bq{z)  are  linear  basic  splines  or  linear  B  splines,  sometimes  referred  to 
as  linear  hat  functions  [Ref.  21].  They  are  defined  as 

for  x  <  Xp-i 


Bp{x)  — 


And: 


B,{z) 


0 

X  -  Xp_i 
Xp  ^Ip— 1 

a:p+i  -  X 

**'p+i  —  ^-p 

0 

0 

Z-Zg-l 

Zq  Zq^i 

Zq+1  -  Z 
Zq+1  Zq 

0 


if  Sp-i  <x<Xp 

iixp  <x  <  Xp^i 
for  X  >  Xp^i 
for  5:  <  z,_i 

if  Zq^i  <  Z  <  Zq 
'A  Zq<Z  <  Zq+I 

for  2:  >  Zp+i 


(51) 


(52) 
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It  is  easier  to  see  how  these  hat  functions  relate  to  each  othe^  when  considering  only 
one  dimension.  Figure  1  shows  that  each  interior  linear  B  spline  only  overlaps  itself 
and  its  two  nearest  neighbors.  This  overlapping  is  an  important  characteristic  which 


will  be  shown  later.  Two  steps  now  need  to  be  done.  First,  define  the  weighting 
functions  shown  in  Equations  (46)  and  (47)  as  Bi(x)  for  the  x  direction,  and  Bj(z) 
for  the  z  direction.  Second,  substitute  the  approximating  (trial)  functions  in  the 
form  shown  in  Equations  (48),  (49),  and  (50)  into  Equation  (42),  and  invoke  the 
Galerkin  procedure.  This  yields, 

0  0 

0  0 

L  K 

+  Jj  T.°°¥oT^'>Bi{x)Bj{z)dzdx 

0  0 

L  K 

=  J  Bi{x)  Bj{z)  dz  dx 

0  0 

L  K 

0  0 
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+  11  Bi{x)Bj{z)  dzdx 

0  0 
L  K 

J  J  Bi{x)  Bj(z)  dz  dx 

0  0 
LK 

Bi{x)  Bj[z)  dz  dx 

0  0 
L  K  ^ 

-  /  /  £  ®  Bi{x)  Bj  (z)  dz  dx 

0  0 
L  A-  „ 

?,(a;)  Bj(z)  dz  dx 

0  0 
L  A-  . 

-J  I  Bi{x)Bj{z)dzdx  (53) 

0  0 

The  streaming  terms  (partial  derivative  terms)  in  Equation  (53)  will  be  integrated 
by  parts.  This  will  accomplish  two  important  results.  First,  since  the  B  splines  are 
only  continuous  in  their  first  derivative,  terms  having  two  partials  operating  on  them 
would  produce  a  value  of  zero.  However,  integrating  by  parts  lowers  the  derivative 
order  applied  to  the  B  splines  by  one.  Thus,  this  legitimizes  their  use.  Second,  but 
probably  the  most  important  result,  is  that  the  integration  procedure  introduces, 
with  relative  ease,  the  natural  boundary  conditions  into  the  finite  elements  [Ref.  13]. 

Integration  by  parts  in  one  dimension  is  defined  in  the  usual  manner, 

b  ^  b 

J  udv  =  uv  —  j  vdu. 


However,  in  two  dimensions,  integration  by  parts  is  done  using  Green’s  theorem: 

JJ  u{V-v)dD  =  I  u{v-n)  dS- JJ  V-  VudD.  (54) 

DSD 

In  three  dimensions,  integration  by  parts  is  defined  as  Gauss’s  theorem  [Ref.  13). 
So  an  area  integral  is  put  in  terms  of  a  surface  integral  using  Green’s  theorem.  To 
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use  Equation  (54),  let, 


V 


A-4.  A 

dx  ^  ^  dz 


k 


n 

V 

u 


1  Tl^  k 

5, -(a:)  Bj{z) 


dz 


n  is. a  unit  vector  normal  to  the  slab’s  surface,  always  poiixting  in  an  outward  direc¬ 
tion.  and  are  the  direction  cosines  of  the  unit  vector  n.  Figure  2  shows  the 
geometry  definitions. 


n  =  k 


n  =  1 


Figure  2:  Two  Dimensional  Slab  Geometry 

If  the  entire  slab’s  surface  is  defined  as  S,  then  performing  the  surface  integra¬ 
tion  in  Equation  (54)  over  the  entire  surface  is  the  same  as  integrating  over  each 
individual  surface  part.  Thus, 

J  dS  =  J  dSi  -b  J  dSi  -b  J  dSz  +  J  dS^. 

S  Si  S2  S3  Si 

The  dot  product  of  V  and  n  then  isolates  the  particular  surface  face  in  ciuestion. 
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Integrating  by  parts  appears  to  complicate  the  equation  because  it  introduces 
more  terms.  However,  later  substitutions  will  be  made  to  reduce  the  number  of 
boundary  terms.  Now  integrating  by  parts  and  rearranging  puts  the  transport  equa¬ 
tion  in  the  form: 


0  0  0  0 


00  ^(N+l) 


B;{x)  Bj{z)  dzdx  =  j  J  Bi{x)  Bj{z)  dz dx 


Jj  J\  U 

+ 1 1  siT^'>  Bi{x)Bj{z)dzdx  + 1 1  v¥^o^Bi{x)Bjiz)dzdx 

0  0  0  0 

L  K  L  K  fiR  (  \ 

+  11  Bi{x)  Bj{z)  dzdx  +  J  J  °  y  Bi{x)  -Qp  dz  dx 
0  0  0  0 

+  J  J  D^^Wo'^^'>Biix)^^Pdzdx  +  l  j  D^^v¥^,^^-^Bj{z)dzdx 


J  J  dzdx  +  jD^° Bi{x)  Bj{z) 


^—11  r  d 


U 

dz-  j  Bi{x)  Bj{z)  ¥^q^  dx 


10  0 


^  K  r  ^ 

D^°Bi{x)  Bjiz)s[^o-^^'>  dx-  I  D^^VBj{z)  Bi{x)^f  dz 

0  0  -0 


JjD^^Bjiz)  Bi{x)s{^i-^^'>  dz 


Equation  (55)  has  seventeen  terms,  including  six  boundary  terms.  The  number  of 
boundary  terms  may  be  reduced  to  two  if  substitutions  derived  from  Equations  (40) 
and  (41)  are  made.  Recalling  that  and  and  making  the 


26 


substitutions  yields, 


0  0  0  0 


-Rii  aB.(x) 


L  K 


0  0 
L  K 


dx  dx 


Bj{z)  dzdx 


+  J  J  Bi(x)  Bj(z)  dzdx  =  j  j  Bi{x)  Bj (z)  dz  dx 


0  0 
L  K 


0  0 


L  K 


+  J  J  si^o-^^'>Bi{x)B j{z)  dzdx  +  :J  J  V  ¥^o^  B{{x)  Bj{z)dzdx 


0  0 
L  K 


0  0 
L  K 


+  11  Bi{x)  Bj(z)  dzdx  +  j  j  Bi{x)  dz  dx 


0  0 
L  K 


0  0 


+ 


J  J  D'Os'S*" Bti=>)^^^dzdx  +  J  J  D"v¥/t>^^^Si(z)dzdx 


0  0 
L  K 


0  0 
K 


^■j  j  D''sf+'>^^Bi(z)dzdx-  j  Bi(z) 


0  0 
L 


S(^+l) 


Bi{x) 


dz 


-  j  Bi{x) 


Bi(z)¥^»^ 


K 


dx 


(66) 


Reducing  the  number  of  boundary  terms  introduces  present  time  step  currents, 
and  which  are  unknown.  However,  we  can  use  Marshak  boundary 

conditions  and  get  these  currents  in  terms  of  9],  [Ref.  11],  [Ref.  15] 

and  [Ref.  16].  For  XZ  geometry,  the  Marshak  condition  for  a  vacuum  is:  [Ref.  11] 


J  n  •  fi Pj”*(cos  0)  cos(m^)  ^(r, E,  ft, t)  dft  =  0  (57) 

nJ2<0 

where  the  angular  flux,  $(r,jF,ft,  f)  is  defined  in  Equation  (22)  and  n  is  a  unit 
normal  vector  at  the  surface  pointing  outwards.  (Equation  (57)  is  for  even  /.) 

Marshak  boundary  conditions  set  the  integral  of  the  incoming  current  to  zero. 
Although  neutrons  that  leave  a  surface  possess  a  finite  probability  of  returning  in 
reality,  the  ones  that  return  usually  have  a  negligible  effect  on  criticality.  It  hcis  been 
determined  that  for  low  order  P„  approximations,  Marshak  boundary  conditions 
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yield  better  results  than  say  Mark  boundary  conditions.  In  fact,  Henry  states  that 
numerical  studies  show  that  Marshak  conditions  lead  to  consistently  more  accurate 
results  up  to  a  Pjg  approximation  [Ref.  1]. 

At  a:  =  0  and  x  =  L,  Bx{x)  =  Ba{x)  =  1.0.  And  at  z  =  0  and  z  =  K, 
Bi{z)  =  Bb{z)  =  1.0.  The  subscripts  A  and  B  on  the  linear  B  splines  indicate 
the  upper  limits  of  the  trial  function  expansions.  If  we  consider  a  Marshak  vacuum 


condition  with  Equation  (57),  then  the  transport  equation  becomes, 

0  0  0  0 
L  K  L  K 

+1 J  B{{x)  Bj{z)  dz  dx  =  J  j  Bi{x)  Bj{z)  dz  dx 


+  /  / 


U 

Bi{x)  Bj{z)  dzdx  +  J  j  V  $^0^  Bi{x)  Bj(z)  dz  dx 


L  K  L  K  Pift  (  \ 

I  J  C^^^Bi{x)  Bj{z)  dzdx  +  J  J  D^°V^[^o'^Bi{x)  ^-^^  dzdx 

00  00 

J  J  D^°S?o-^^^Bi{x)^^^dzdx  +  J  J  D^^vMf^^^Bjiz)dzdx 
0  0  0  0 
L/  K  K 

/  /  iJ"  sS"+'>  Biiz)  dzdx-l  jf:  B,(z)  Biiz)  dz 

n  n  n  <7=1 


1 


-\jt,  B,(z)  BMdz  -  5  /  E  BpW  Bi{x)  dx 

0  0 

-  \  j  BA^)Bi(x)dx 

n  P=t 


Recall  that  the  expansion  of  the  flux  in  the  x  direction  went  from  p  =  1  to  p  = 
A.  In  the  z  direction,  the  expansion  went  from  q  =  1  to  q  =  B.  Therefore,  in 
Equation  (58),  the  boundary  flux  terms  are  defined  as,  (see  Figure  2  on  page  25) 


*00l7 

is  on  the  surface  where 

X  =  0 

z  =  0  to  K 

*004, 

is  on  the  surface  where 

X  =  L 

z  =  0  to  K 

*00pi 

is  on  the  surface  where 

z  =  Q 

X  =  0  to  L 

*00pij 

is  on  the  surface  where 

II 

X  =  0  to  L 

Now  the  whole  transport  equation  will  be  expanded: 


0  0  P=1  9=1 

n  n  P—1  9=1 


0  0  P-^  9= 


L  K 


A  B 


+  /  /  S°°  E  E  ^^mt^B,(x)B,(z)B,(x)Bi(z)dzdx 

0  0  p~^  9”^ 

L  K  A  B 

=  /  j  ^^EE^To*'^B,{x)  B,[z)B,(x)Bj{z)dzdx 

0  0  P~^  9=1 


L  K 


A  B 


+ 


J  J  Bp{x)  Bg{z)  Bi{x)  Bj{z)  dz  dx 

0  0  P~^  9=1 


L  K 


A  B 


+  /  J  ^  fill  '^i^olgBp{x)Bg{z)Bi{x)Bj{z)dzdx 

0  0  9=1 

L  K  J9 

+  //  llflC^^'^Bp{x)Bg{z)Bi{x)Bj{z)dzdx 

0  0  9=1 

+  /  /  D^°y  E  E  ^'■;2,,B,(x)B,(z)Bi(x)^^dzdx 

5  5  P=1  7=1 

+  j  Ib'°  EE  S%t''B,[x)B,(z)Bi(x)^^dzdx 

0  0  p=l  y— 1 

+  jjD"VEE  B,{x)  B,{z)  ^  W  <i.  dx 


p=l  9=1 
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0  0  P=^  9=1 

-5/E*io^’-B,W®iW* 

^  0  9=1 

^  0  9=1 
1  ■1'  A 

^  0  P=1 

-  5  /  E  *0?;,'’  BpW  5i(^)  (59) 

As  previously  stated  and  shown  in  Figure  1,  each  linear  5  spline  overlaps  itself  and 
its  two  nearest  neighbors.  This  implies  that  p  =  i  —  +  l  and  q  =  j  —  l,i,i  +  !• 

Also,  the  integration  may  now  be  over  each  element,  since  the  boundary  fluxes  are 
set,  and  each  element  carries  the  boundary  terms  with  them  [Ref.  13].  e*  will  denote 
integration  of  the  elements  in  the  x  direction  and  will  likewise  denote  integration 
of  the  elements  in  the  ^  direction.  Therefore, 


11°'°%  E  ^fC^B,{x)^-^Bix)^-^dzdx 

ex  cz  p=t-l  q=j~l 

+  11^"%  E  C;'^B,{x)^Biix)dzdx 

cx  Cz  p=«-l  g=j— 1 

r  f  i+1  i+1 

+  /  /s“”  E  E  ^?o*;^B,(x)B,(.z)BAx)Si{z)dzdx 

Cx  Cx  P=t-1  g=i— 1 

.  -  t+1  i+1 


Cx  Cz 


m  m  ^  r* 

+  /  /  E  E  Si'i*»B,{x)B,iz)B,{x)Bi{z)dzdx 
ix  iz  p=«'-i  9=i-i 
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t+1  i+1 


r  r  *  *  * 

+  /  /  E  E 

Cx  c-  P=*— 1  9=i-l 

+  /  j  '■£  ’£  9\'ilB,(x)B,(z)Bi{x)^^dzdT 

Cx  ix  P=i-1  g=j-l 

Cl  Cl  P=i— 1  q=i— 1  ^ 

+  ff  D''V  E  £ 

Cx  iz  P=«— 1  q=}~\  ^ 


Cx  c-  P=*-l  9=i-l 

i+1 


ax 

i+1 


1  ^  jTl  1  .  J-t-i 

-5/  E  ^ml''B,{x)Bi{z)dz---f  E 

e*  9=i-l  ^  Cl  9=i-l 

-5/  E  E 

Cx  P=‘-l  '‘e,  P=«“l 


(60) 


The  integration  summaries  for  the  elements  are  given  in  Appendix  B.  Each  mate¬ 
rial  matrix  is  assumed  to  be  piece- wise  continuous  in  each  xz  interval.  Bringing  the 
boundary  terms  over  to  the  right  side  of  the  equation  and  integrating  Equation  (60) 
with  Kronecker  delta  notation  yields. 
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3.6  Coefficient  Matrix  Example 

When  Equation  (61)  is  implemented  a  square  symmetric  positive  definite  block  ma¬ 
trix  equation  is  produced.  The  matrix  equation  is  in  the  form  of  =  S. 

A  is  the  coefficient  matrix,  are  the  expansion  coeflScients  that  physically 

represent  the  total  flux  at  mesh  points  ij,  and  S  is  the  source  vector.  If  we  are 
considering  a  system  that  has  G  energy  groups,  then  each  entry  in  A,  A{j,  is  a  GxG 
matrix.  Likewise,  each  of  the  entries  in  and  S  are  Gxl  vectors.  Thus 

we  have  a  block  matrix  equation.  Block  multiplication  can  be  carried  out  in  the 
same  way  as  traditional  matrix  multiplication  since  all  entries  will  have  the  proper 
dimensions  [Ref.  20]. 

Suppose  we  have  a  slab  with  16  mesh  points,  4  in  the  x  direction  and  4  in  the 
2  direction  as  shown  in  Figure  3.  By  applying  Equation  (61)  for  the  sixteen  mesh 
points,  we  will  get  a  16  by  16  matrix,  A,  and  16  by  1  vectors,  and  S.  This 

system  is  shown  in  Figure  4.  The  block  matrix  has  a  bandwidth  of  11  with  nine 
nonzero  diagonals.  For  this  example,  there  are  100  nonzero  A,j  entries.  Each  entry 
is  a  material  defined  matrix  for  an  interval  along  with  factors  resulting  from 
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the  finite  element  integrations.  For  a  particular  spatial  interval,  the  material  cross 
sections  are  assumed  to  be  piece-wise  constant. 

This  development  of  the  coefficient  matrix  was  done  in  a  “brute  force”  way. 
Traditional  finite  element  schemes  assemble  the  coefficient  matrix  in  a  much  simpler 
way  by  using  a  local  to  global  mesh  point  numbering  scheme.  This  entails  the 
addition  of  local  coefficient  matrices  to  yield  the  global  coefficient  matrix.  The 
source  vector  is  likewise  developed.  Thus  one  advantage  of  finite  element  schemes 
was  not  taken  advantage  of  here.  But,  it  is  clearer  to  see  how  each  entry  in  the 
coefficient  matrix  was  derived  by  developing  the  coefficient  matrix  in  this  manner. 

Now  we  summarize  this  process  for  RZ  geometry.  Basically,  everything  is  done 
the  same  as  for  XZ  geometry,  but  the  RZ  Ccise  is  slightly  more  complicated  because 
of  its  streaming  terms. 
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Figure  4:  Spatial  Matrix  Example 


4  Solution  in  RZ  Geometry 


We  now  do  much  the  same  procedure  to  obtain  an  approximate  solution  to  the 
transport  equation  in  RZ  geometry.  The  cylindrical  coordinate  system  depends  on 
azimuthal  angle,  Xj  a  distance,  r,  measured  from  the  z  axis,  and  the  coordinate  z 
mecisured  on  the  z  axis  itself  from  the  x  y  plane.  Figure  5  shows  the  geometry  defini¬ 
tions.  If  we  have  azimuthal  symmetry,  then  the  neutron  density  is  not  changing  with 
respect  to  the  azimuthal  angular  coordinate  x-  Therefore  the  rate  of  change  of  the 
neutrons  along  a  streaming  path  is  accurately  described  by  Equation  (8).  The  flux 
and  source  may  be  expanded  as  before  in  Equations  (2)  and  (3).  Thus  implementing 


Figure  5:  Cylindrical  Geometry 

Equations  (2),  (3),  and  (8)  into  Equation  (1)  puts  the  transport  equation  in  a  form 
much  like  Equation  (12): 
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(62) 


The  angular  dependency  of  the  angular  neutron  flux  is  described  the  same  as  for  XZ 
geometry.  Spherical  harmonics  have  been  used  again  for  the  angular  discretization 
of  the  direction  vector  Cl  because  they  completely  describe  the  angular  dependency 
of  Cl.  Note,  because  of  the  relationship,  P;“’"(cos//)  =  (— I)*"  P,’"(cos /?),  it 

is  unnecessary  to  include  negative  m  [Ref.  22]  [Ref.  23].  To  proceed,  we  need  two 
more  equations  in  addition  to  Equation  (16).  They  are. 


^  cos(m^)  =  —m  sin(m^) 

1  1 

s\n{m(l>)  sin  ^  =  -  cos(m  —  1)^  —  -  cos(m  +  l)<f> 


(63) 

(64) 


As  before  we  assume  fission  to  be  an  isotropic  event,  invoke  the  addition  theorem, 
and  rearrange.  The  transport  equation  now  takes  the  form: 
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(65) 


As  with  XZ  geometry,  Equation  (65)  will  be  manipulated  by  implementing  the  re¬ 
currence  relationships,  adjusting  the  indices  to  put  it  back  in  terms  of  spherical 
harmonics,  and  then  multiplying  by  Pj^ (cosO)  cos(N^)  and  integrating  in  all  direc¬ 
tions.  With  stipulated  values  for  K  and  N  and  using  Euler’s  backward  differencing 
scheme  for  time  discretizations,  the  following  vector  equations  are  produced: 
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(66) 
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or 

Now  we  substitute  the  value  of  from  Equation  (67)  and  the  value  of 

from  Equation  (68)  into  Equation  (66)  to  yield, 
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Equation  (69)  is  the  RZ  equivalent  to  Equation  (42)  for  XZ  geometry,  but  Equa¬ 
tion  (69)  has  more  terms,  and  some  terms  have  a  ^  factor  in  them.  As  the  finite 
element  spatial  discretization  is  applied  to  Equation  (69),  some  of  the  extra  terms 
cancel  with  terms  produced  when  integrating  by  parts.  We  now  proceed  with  the 
finite  element  discretization. 


4.1  RZ  Finite  Element  Discretization 


The  neutron  sources  in  Equation  (69)  may  be  expanded  with  the  same  basis  hat 
functions  in  the  r  and  2  directions  as  was  done  in  the  XZ  geometry  in  Equations  (48), 
(49),  and  (50).  The  hat  functions  in  the  z  direction  are  defined  as  in  Equation  (52). 
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The  functions  in  the  r  direction  are  defined  as: 


Bp{r) 
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if  J’p-i  <r  <rp 
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rp+i  -  rp 

0  for  r  >  Tp+i 
Now  substituting  the  trial  functions  into  Equation  (69)  yields: 
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(71) 


The  Galerkin  procedure  can  now  be  implemented.  However,  integration  over  the 

domain  is  now  defined  for  a  cylinder  with  a  radius  from  0  to  il  and  a  height  from  0 

K  R 

to  K.  Therefore,  the  domain  integration  is,  /  /  2zrdrdz.  Implementing  this,  the 

0  0 

transport  equation  is  now: 

K  R 
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0  0 

Green’s  theorem  is  again  used  to  integrate  the  streaming  terms  by  parts.  This  re¬ 
duces  the  continuity  requirement  of  the  trial  functions  and  implements  the  natural 
boundary  conditions.  And  since  each  element  carries  the  boundary  conditions  with 
them,  the  integrations  may  be  designated  over  each  element.  Let  indicate  inte¬ 
gration  over  the  elements  in  the  r  direction  and  let  Cz  indicate  integration  over  the 
elements  of  the  2  direction.  Implementing  the  Galerkin  procedure,  the  transport 


equation  is  now: 
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Appendix  B  gives  the  integrating  summaries.  Again,  Marshak  v<acuum  boundary 
conditions  are  implemented  to  substitute  for  the  and  coefficients  in 
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terms  of  coefficients.  Using  the  Kronecker  delta  where  p  =  i  —  l,z,  and  z  + 1 

and  q  =  j  —  1,3 1  and  j  + 1  produces: 
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This  equation  will  also  produce  a  block  matrix  equation  with  nine  diagonals.  Now 


we  will  address  the  solutio’  *^echniques  to  solve  Equations  (61)  and  (74). 
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5  Solution  Technique 

Using  Equation  (61)  or  (74)  produces  a  block  matrix  equation  as  show  in  the  example 
in  Figure  4.  The  coefficient  matrix,  A,  will  always  have  nine  full  diagonals.  For  the 
16  by  16  example  shown,  there  were  100  nonzero  entries  in  A.  That  means  that 
156  entries  contain  nothing  but  zeros.  Recall  that  each  entry  in  A  is  itself  a  G  by  G 
matrix  representing  G  energy  groups.  Storage  of  all  those  zero  matrices  would  be 
very  inefficient.  Therefore,  FMP2DT  uses  some  storage  schemes  to  avoid  using 
excessive  computer  memory.  And  because  this  system  can  be  extremely  sparse  and 
large,  depending  on  the  number  of  mesh  points  and  energy  groups  selected,  FMP2DT 
used  some  special  algorithms  to  'calculate  a  solution.  The  following  is  information, 
including  some  background,  concerning  FMP2DT’s  solution  algorithm. 

If  A  G  R”’'”  is  a  symmetric  positive  definite  matrix,  then  a  lower  triangular 
matrix  G  G  R”*”  exist  with  positive  diagonal  entries  such  that  A  =  GG^.  Splitting 
A  like  this  is  known  as  the  Cholesky  decomposition  of  the  matrix  A.  To  solve  the 
system  Ax  =  b  using  a  Cholesky  algorithm  entails  computing  A  =  GG^  and  then 
solving  Gy  =  b  and  then  G^x  =  y  [Ref.  18].  This  solution  technique  is  stable  and 
efficient  for  solAng  large  banded  systems.  This  solution  technique  is  called  factoring 
the  coefficient  matrix  A,  and  is  referred  to  as  a  direct  method  of  solution. 

However,  for  very  large  and  sparse  systems,  which  could  be  our  case,  direct 
methods  are  often  not  efficient  enough.  For  a  linear  system,  iterative  methods  are 
more  suitable.  One  iterative  method  is  call  the  conjugate  gradient  method  [Ref.  18). 
It  involves  minimizing  a  functional  <f>{x)  such  that, 

=  ix^  A  X  —  x^  b. 

If  A  is  symmetric  and  positive  definite,  then  minimizing  the  above  expression  is 
the  same  as  solving  for  Ax  =  b  [Ref.  18].  But  convergence  of  a  steepest  descent 


algorithm  may  be  extremely  slow.  Therefore,  preconditioning  A  is  desirable  to  speed 
convergence.  One  preconditioning  strategy  is  developing  an  incomplete  Cholesky 
factorization  of  A.  This  involves  calculation  of  some  lower  triangular  matrix  that  is 
somewhat  close  to  the  actual  Cholesky  lower  triangle  matrix  G  [Ref.  18]. 

This  solution  technique  FMP2DT  uses  to  do  this  is  the  incomplete  Cholesky 
conjugate  gradient  algorithm  [Ref.  33].  Two  iteration  schemes  are  used.  The  inner 
iterations  are  over  the  spatial  mesh  points.  The  outer  iterations  are  over  the  specified 
energy  groups. 

Now  we  examine  one  of  the  important  features  of  FMP2DT.  FMP2DT  has  the 
ability  to  calculate  a  delayed  x  spectrum  if  that  spectrum  is  not  known.  If  it  is 
known,  then  it  may  be  entered  in  the  input  deck.  However,  in  most  cases,  these  val¬ 
ues  are  not  easily  found.  Therefore,  the  next  section  shows  how  FMP2DT  calculates 
these  values  for  the  user. 
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6  Calculation  of  %f 


In  the  multigroup  derivations,  Equation  (27)  has  x  parameters  for  both  prompt 
and  delayed  neutrons.  Most  cross  section  sets  contain  information  defining  Xp:  but 
they  do  not  define  Xk-  values  for  xi  are  known,  then  they  may  be  entered  in 
FMP2DT’s  15iWr  array  in  the  input  deck.  If  this  array  is  filled  with  zeros,  and 
FMP2DT  is  running  a  fission  problem,  then  an  internal  flag  is  triggered  that  causes 
FMP2DT  to  assume  that  the  user  does  not  know  these  values.  FMP2DT  then 
proceeds  to  calculate  them  and  continues  the  solution  process. 

The  spectra  for  prompt  neutrons  are  easily  found  in  literature.  There  are  many 
formulae  that  have  been  derived  to  fit  the  data  mathematically  by  Watt,  Cranberg, 
and  others  [Ref.  34]  [Ref.  35).  Each  analytical  fit  is  due  to  examination  of  data 
over  certain  energy  ranges  [Ref.  36].  FMP2DT  assumes  a  Maxwellian  distribution 
function  that  is  defined  as  [Ref.  30], 


kT  IkT 
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-B 

e*i- 


(76) 


where  k  is  Boltzmann’s  constant  and  T  is  the  spectrum  temperature  (K).  The  av¬ 
erage  energy  is  given  by 

E  =  \kT. 

2 

However,  the  delayed  neutron  spectrum  is  not  as  well  established.  Several  dif¬ 
ferent  data  exist  with  various  uncertainty.  FMP2DT  assumes  that  delayed  neu¬ 
trons  follow  the  same  sort  of  Maxwellian  distribution  as  prompt  neutrons  [Ref.  30]. 
Assuming  that  kT=1.29  for  prompt  neutrons,  or  that  kT=0.29  for  delayed  neu¬ 
trons  [Ref.  31],  Equation  (75)  may  be  used  to  calculate  either  a  prompt  or  delayed 
spectrum.  (The  values  for  T  above  are  for  and  kT  is  in  MeV.) 


For  a  X  spectrum,  we  note  that 

CO 

J  x{E)dE=  1.0. 

0 

For  a  multigroup  approximation,  x^  is  calculated  as 

Eg  — I 

X‘=  j  x(E)dE. 

Eg 

FMP2DT  calculates  Xk  likewise  using  an  adaptive  quadrature  scheme  based  on 
Gauss-Kronrod  algorithms.  And  as  stated,  changing  the  value  for  kT  can  cause 
an  evaluation  for  a  prompt  x  calculation.  Also,  if  a  better  distribution  function 
is  desired,  then  a  short  function  subprogram  may  be  added  to  the  code  and  the 
quadrature  scheme  can  evaluate  it. 

Using  a  47  neutron  group  structure  shown  in  Tables  1  and  2,  the  following 
results  were  generated.  First,  a  value  of  ^'r=1.29  was  used  to  generate  a  prompt  x 
distribution  shown  in  Figure  6.  Note  that  although  the  neutron  groups  go  to  17.33 
MeV,  the  graph  stops  at  10  MeV.  This  is  because  the  Xp  s  are  approximately  zero 
beyond  that  point.  Using  FMP2DT’s  quadrature  scheme  with  a  value  of  kT=0.29, 
Figure  7  was  generated  for  the  delayed  x^-  It  was  truncated  at  only  2.231  MeV 
because  of  the  same  reason  as  the  prompt  data. 

It  is  important  to  notice  that  the  delayed  x  s  peak  at  lower  energy  groups  than 
the  prompt  x  s  do.  Also  the  prompt  x  s  range  over  a  much  larger  number  of  groups. 
This  is  expected  since  prompt  fission  neutrons  tend  to  be  born  with  higher  energies 
than  the  delayed  fission  neutrons. 

Now  we  want  to  establish  FMP2DT’s  computational  integrity.  First,  a  compar¬ 
ison  between  flux  shapes  was  made  with  a  flux  calculation  by  a  two  dimensional, 
two-group,  space-time  diffusion  code  called  TWIGL.  Then  FMP2DT  was  bench- 
marked  with  some  exact  flux  calculations  in  both  slab  and  cylindrical  geometries. 
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:  Group 

Upper  Bound  (MeV) 

Lower  Bound  (MeV) 

1 

1.7330E+01 

1.4190E+01 

2 

1.4190E+01 

1.2210E+01 

3 

L2210E+01 

l.OOOOE+01 

4 

l.OOOOE+01 

8.6070E+00 

5 

8.6070E+00 

7.4080E+00 

6 

7.4080E+00 

6.0650E+00 

7 

6.0650E+00 

4.9650E+00 

8 

4.9650E+00 

3.6780E+00 

9 

3.6780E+00 

3.0110E+00 

10 

3.0110E+00 

2.7250E+00 

11 

2.7250E+00 

2.4660E+00 

12 

2.4660E+00 

2.3650E+00 

13 

2.3650E+00 

2.3450E+00 

14 

2.3450E+00 

2.2310E+00 

15 

2.2310E+00 

L9200E+00 

16 

1.9200E+00 

L6530E+00 

17 

1.6530E+00 

1.3530E+00 

18 

1.3530E+00 

1.0020E+00 

19 

1.0020E+00 

8.2080E-01 

20 

8.2080E-01 

7.4270E-01 

21 

7.4270E-01 

6.0810E-01 

22 

6.0810E-01 

4.9780E-01 

23 

4.9780E-01 

3.6880E-01 

24 

3.6880E-01 

2.9720E-01 

Table  1:  Part  1:  47  Neutron  Group  Structure 


54 


Group 

Upper  Bound  (MeV) 

Lower  Bound  (MeV) 

25 

2.9720E-01 

1.8310E-01 

26 

1.8310E-01 

l.llOOE-01 

27 

l.llOOE-01 

6.7370E-02 

28 

6.7370E-02 

4.0860E-02 

29 

4.0860E-02 

3.1820E-02 

30 

3.1820E-02 

2.6050E-02 

31 

2.6050E-02 

2.4170E-02 

32 

2.4170E-02 

2.1870E-02 

33 

2.1870E-02 

1.5030E-02 

34 

1.5030E-02 

7.1010E-03 

35 

7.1010E-03 

3.3540E-03 

36 

3.3540E-03 

1.5840E-03 

37 

L5840E-03 

4.5400E-04 

38 

4.5400E-04 

2.1440E-04 

39 

2.1440E-04 

1.0130E-04 

40 

1.0130E-04 

3.7260E-05 

41 

3.7260E-05 

1.0670E-05 

42 

1.0670E-05 

5.0430E-06 

43 

5.0430E-06 

1.8550E-06 

44 

1.8550E-06 

8.7640E-07 

45 

8.7640E-07 

4.1390E-07 

46 

4.1390E-07 

9.9990E-08 

47 

9.9990E-08 

l.OOOOE-11 

Table  2:  Part  2:  47  Neutron  Group  Structure 
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Figure  6:  Prompt  Group  x  for 
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Figure  7:  Delay  Group  x  for 
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7  Code  Benchmark 


While  there  are  many  publications  of  time  dependent  neutron  transport  work  avail¬ 
able,  it  is  extremely  difficult  to  find  a  two-dimensional  problem  with  a  delayed 
neutron  source  so  that  FMP2DT  can  be  benchmarked.  While  problems  of  this  na¬ 
ture  are  often  cited,  the  published  results  usually  do  not  contain  enough  information 
about  parameters  and  cross  sectional  data  to  reproduce  them.  In  1968,  a  technical 
report  was  published  showing  results  for  a  time-dependent,  two  dimensional  slab 
problem  solved  with  a  code  called  TWIGL  [Ref.  28].  TWIGL  is  a  two  dimensional, 
two-group,  space-time  diffusion  equation  solver  that  incorporates  temperature  feed¬ 
back. 

TWIGL  Wcis  used  to  compare  FMP2DT’s  fast  flux  shape.  To  benchmark  FMP2DT’s 
computation  accuracy,  data  in  a  report  by  B.  D.  Ganapol  [Ref.  29]  was  used.  First 
we  show  the  results  for  the  TWIGL  comparison,  and  then  two  benchmarks  using 
Ganapol’s  data  for  infinite  RZ  and  XZ  geometries. 

7.1  TWIGL  Comparison 

TWIGL  was  used  to  compare  FMP2DT’s  fast  flux  shape.  This  offered  at  least 
some  sort  of  comparison  between  source  vectors  for  the  two  codes,  and  both  source 
vectors  have  a  delayed  neutron  contribution  in  them.  TWIGL  is  a  diffusion  code. 

It  basically  solved  the  following  equations: 

V  •  Di (r,  i)  V  <f)i (r,  t)  -  S, (r,  1)  t) 

+  {l-P)  [j/S/,(r,i)  -f  i/S/j(r,0Mr,t)] 

i=i 


(76) 


1 

V-Z)2(r,i)  V^2(r,f)  -S2(r,f)^2(r,i)  +  Sr,(r,f)^i(r,f)  -  —  (77) 

V2  eft 

and 

Q 

^  C((r,  <)  =  [j'Sy,  (r,  i)  (r,  i)  +  i-B/,  (r,  i)  ^2(r,  i)]  -  Af Oi(r,  i),  i  =  1,  ■  •  ■  / 

(78) 

where  r  represents  x,^  for  slab  geometry.  The  slab  geometry  for  this  problem  is 
shown  in  Figure  8.  The  equations  above  are  solved  by  TWIGL  subject  to  zero  flux 
boundary  conditions  on  all  external  surfaces.  At  time  t  =  0,  the  reactor  is  critical,  i.e. 
ken  =  1-0.  The  initial  flux  is  calculated  using  this  steady  state  condition.  TWIGL 


Figure  8:  TWIGL  Slab  Geometry 

discretizes  the  time  dependent  flux  using  a  backward  differencing  scheme  (for  this 
problem)  and  a  central  differencing  scheme  for  the  precursor  terms. 

For  most  problems,  FMP2DT  reads  cross  sectional  data  from  an  input  tape. 
However,  data  can  be  input  directly.  The  input  parameters  used  were  derived  from 
the  data  in  the  TWIGL  report.  Since  TWIGL  is  a  diffusion  code,  the  values  cor¬ 
responding  to  the  entries  in  FMP2DT’s  input  deck  are  set  to  zero.  This  is 
known  as  the  diffusion  approximation  for  the  Pi  calculations.  This  is  also  done  for 


the  calculation  of  the  initial  flux  at  time  t  =  0,  which  wais  calculated  by  FEMP2D. 
FMP2DT  can  calculate  its  own  initial  flux  by  using  a  very  large  At.  This  will  cause 
all  the  terms  to  be  approximately  zero.  From  Appendix  A,  it  is  seen  that 
all  the  time  dependent  terms  will  drop  out.  But  using  an  initial  flux  calculated 
by  FEMP2D  is  more  cost  effective  since  a  large  At  in  FMP2DT  still  entails  some 
unnecessary  calculations  .  like  V  =  0. 

Table  3  compares  the  TWIGL  and  FMP2DT  initial  fluxes.  Because  of  symmetry, 


X(cm) 

TWIGL 

FEMP2D 

10 

6.26847E+13 

5.9368E+13 

20 

1.92851E+14 

1.8581E+14 

30 

5.34671E+14 

5.5854E+14 

40 

9.37259E+14 

9.4485E+14 

50 

1.08474E+15 

1.0848E+15 

60 

9.39969E+14 

9.4610E+14 

70 

5.39333E+14 

5.6072E+14 

80 

2.00865E+14 

1.9003E+14 

90 

8.28086E+13 

7.1108E+13 

100 

5.38778E+13 

4.3079E+13 

110 

8.28086E+13 

7.1108E+13 

120 

2.00865E+14 

1.9003E+14 

130 

5.39333E+14 

5.6072E+14 

140 

9.39969E+14 

9.4610E+14 

150 

1.08474E+15 

1.0848E+15 

160 

9.37259E+14 

9.4485E+14 

170 

5.34671  E+14 

5.5854E+14 

180 

1.92S51E+14 

1.8581E+14 

190 

6.26847E+13 

5.9368E+13 

Table  3:  TWIGL  and  FMP2DT  Initial  Fast  Fluxes 


the  flux  data  are  for  z  ^  14.142  cm.  There  are  no  flux  values  given  in  Table  3  for 
X  =  0  or  X  =  200  cm  since  TWIGL  sets  these  values  to  zero.  FMP2DT,  however, 
does  not  set  these  fluxes  to  zero  because  it  models  a  vacuum  boundary  condition. 
So  the  two  codes  should  agreo  more  for  the  interior  mesh  point  flux  calculations. 
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Figure  9  shows  the  plot  of  the  TWIGL  flux  at  time  i  =  0. 

The  material  mean  free  paths  (MFP)  are  calculated  in  Appendix  C  and  displayed 
in  Table  4.  It  is  obvious  that  mesh  spacing  with  Az  =  14.142  cm  and  Ax  =  10  cm, 
which  was  used  in  the  TWIGL  report,  is  much  larger  than  these  MFPs.  This  suggest 
that  the  thermal  flux  calculations  could  be  suspect.  Indeed,  both  the  FEMP2D  and 
FMP2DT  calculations  showed  that  the  thermal  flux  shape  varied  greatly  as  finer 
mesh  spacing  was  chosen.  The  fast  energy  group,  with  its  longer  MFPs,  was  the 
least  affected  by  mesh  spacing.  Figure  10  shows  the  initial  FMP2DT  flux  calculated 
by  FEMP2D.  It  was  done  using  the  reported  TWIGL  mesh.  To  be  sure  that  this 
mesh  spacing  was  sufficient  to  define  the  fast  flux,  another  FEMP2D  run  was  made 
using  twice  the  reported  TWIGL  mesh  points.  (The  TWIGL  Aa:  and  Az  values 
where  cut  in  half.)  Figure  11  shows  this  result.  There  is  no  substantial  cliange 
between  the  two  FEMP2D  fluxes. 


Group 

Materials  1  &3 

Material  2 

1 

4.1701  cm 

3.5702  cm 

2 

1.5601  cm 

2.1000  cm 

Table  4:  Material  Mean  Free  Paths 


TWIGL  set  the  initial  condition  of  the  slab  to  be  critical.  Using  TWIGL  mesh 
spacing,  FEMP2D  calculated  a  =  1.01  IIS,  and  using  twice  as  many  mesh  points, 
FEMP2D  calculated  a  kea  =  1.026S5S.  So  the  initial  calculations  for  the  slab  arc 
very  close  for  both  codes.  Therefore,  the  TWIGL  mesh  seems  to  be  good  enough  to 
define  the  fast  flux. 
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7.1.1  Transient  Calculations 

The  material  cross  sections  changed  with  time  due  to  two  causes.  The  first  has 
to  do  with  a  linearly  changing  cross  section.  Material  1  in  Figure  8  differs  from 
material  3  because  material  1  has  a  time  dependent  thermal  absorption  shown  in 
Figure  12.  The  second  is  due  to  temperature  change.  The  TWIGL  code  assumes  a 


Figure  12:  Time  Dependent  Thermal  Absorption 

coolant  flow  along  the  z  axis.  Since  it  sets  the  flux  to  zero  on  the  slab  surface,  the 
coolant  has  no  direct  neutronic  elect  such  as  absorption,  reflection,  etc.  However,  it 
couples  with  a  fission  power  calculation  to  establish  a  core  and  coolant  temperature. 
The  material  cross  sections  are  then  adjusted  for  the  change  in  temperature  after 
each  time  step  convergence.  The  time  interval  which  the  TWIGL  had  the  smallest 
temperature  change  was  chosen  for  this  comparison.  For  time  t  =  0  to  t  =  0.01  sec 
there  was  no  change  in  the  core  temperature.  However,  there  was  a  small  change  in 
the  coolant  temperature.  This  is  shown  in  Table  5.  FMP2DT  could  have  modeled 
this  change  if  the  TWIGL  report  had  given  the  temperatures  at  the  end  of  each 
time  step.  However,  it  did  not.  But  the  error  should  be  very  small  since  the  coolant 
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Zone  1 

Zone  2 

Zone  3 

Zone  4 

Zone  5 

0.00 

500.5680  °F 

513.3180  "F 

500.6400  °F 

513.3180  °F 

500.5680  °F 

0.01 

500.5681  °F 

513.3188  °F 

500.6399  °F 

513.3181  °F 

500.5680  °F 

Table  5:  Coolant  Temperatures  by  Zone 


temperature  change  for  each  zone  is  not  significant. 

The  time  dependent  thermal  absorption  for  material  1  can  be  d^cribed  mathe¬ 
matically  as, 


m = 


0.43  For  t  <  0.005  sec 
0.44  —  2t  For  0.005  <t  <  0.01  sec 


(79) 


Except  for  a  small  change  due  to  temperature,  the  D2  TWIGL  parameter  re¬ 
mains  constant  in  the  interval.  Therefore,  Sj  for  the  FMP2DT  calculations  remains 
constant.  Since  Sj  =  -j-  must  be  time  dependent.  This  can  be 

expressed  mathematically 


y7-*2(f  \  _  I  0.211 

*0  1  0.201-h2i 


For  t  <  0.005  sec 
jL-or  0.005  <t  <  0.01  sec 


(80) 


Using  Equations  (79)  and  (80)  above  the  cross  section  value  for  each  parameter  can 
be  calculated  for  each  time  step.  Then  the  value  calculated  for  the  present  time 
step  must  be  subtracted  from  the  value  for  that  parameter  used  during  the  Icist  time 
step  to  obtain  some  AS.  That  AS  goes  into  the  input  deck  in  the  20**,  21**,  and 
22**  arrays  (shown  in  Appendix  C).  If  a  cross  section  has  no  time  dependency,  then 
its  entries  will  be  zero.  (Note,  the  initial  input  into  these  arrays  must  contain  the 
steady  state  values.) 

The  TWIGL  flux  at  time  t  =  0.01  sec  is  given  in  Table  6.  The  peaks  occur  at 
a;  =  50  cm.  and  x  =  150  cm,  or  in  material  zones  2  and  4  which  are  composed  of 
material  1  and  material  3  respectively.  Figure  13  shows  the  TWIGL  flux  and  the 
two  distinct  peaks. 

The  two  peaks  differ  because  of  the  change  in  the  thermal  absorption  in  ma- 
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X(cm) 

10 

<!> 

1.835550-1-14 

20 

5.668937-1-14 

30 

L578674-fl5 

40 

2.774679-1-15 

50 

3.211949-fl5 

60 

2.779915-M5 

70 

1.587650-fl5 

80 

5.824532-H4 

90 

2.22951 1-f  14 

100 

1.039940-fl4 

no 

1.038213-fl4 

120 

2.137265-1-14 

130 

5.582470-1-14 

140 

9.682363-1-14 

150 

1.115436-1-15 

160 

9.629778-1-14 

170 

5.491288+14 

180 

1.980289+14 

190 

6.435757+13 

Table  6:  TWIGL  Flux  at  Time  t  =  0.01  sec 

terial  1.  Thermal  neutrons  have  a  larger  probability  of  inducing  fission  than  fast 
neutrons.  Neutrons  emitted  as  a  result  of  a  fission  event  are  high  energy,  or  fast 
neutrons.  Therefore,  material  absorption  of  thermal  neutrons  result  in  a  decline  of 
fission  events.  Conversely,  less  absorption  increases  the  thermal  neutron  popula¬ 
tion,  and  increases  fission  events.  More  fission  events  then  increase  the  fast  neutron 
population. 

It  follows  then  that  since  the  thermal  absorption  cross  section  decreases  in  ma¬ 
terial  1,  more  fission  events  occur  there.  This  increases  the  fast  neutron  population, 
anf  causes  a  larger  peak  in  material  1. 

FMP2DT  should  also  display  the  same  shape.  Figure  14  shows  this  to  be  the 
case.  Figure  15  shows  the  results  using  twice  the  reported  TWIGL  mesh  spacing. 


In  both  cases,  FMP2DT’s  flux  shape  was  very  similar  to  TWIGL’s. 

FMP2DT  made  twelve  runs  to  insure  that  its  flux  shape  was  correctly  defined. 
Table  7  shows  the  different  configurations  modeled.  Runs  1  through  4  were  done 
using  the  reported  TWIGL  mesh  spacing.  Note,  that  for  runs  1  and  2  there  are 
differences  in  the  number  of  time  intervals  and  At  s.  The  parameters  for  run  1  are  the 
same  as  those  that  TWIGL  used  for  its  calculations.  TWIGL  used  one  time  interval 
with  ten  time  steps  to  model  from  t  =  0  to  t  =  0.01  sec,  with  each  At  =  0.001  sec. 
From  Figure  12  it  is  shown  that  between  times  t  =  0  and  t  =  0.005  sec,  there  is  no 
change  in  the  thermal  absorption  cross  section.  TWIGL  still  used  five  time  steps 
there  even  though  no  physical  process  was  changing. 


Run 

Intervals 

Ati  (sec) 

At2  (sec) 

Ax(cm) 

A2(cm) 

Total  Steps 

1 

1 

0.001 

N/A 

10.0 

14.142 

10 

2 

2 

0.005 

0.001 

■CQI 

14.142 

6 

3 

2 

0.005 

0.0005 

10.0 

14.142 

11 

4 

2 

0.005 

0.00025 

10.0 

14.142 

21 

5 

1 

0.001 

N/A 

20.0 

14.142 

10 

6 

2 

0.005 

0.001 

20.0 

14.142 

6 

7 

2 

0.005 

0.0005 

20.0 

14.142 

11 

8 

2 

0.005 

0.00025 

20.0 

14.142 

21 

9 

1 

0.001 

N/A 

5.0 

7.071 

10 

10 

2 

0.005 

0.001 

5.0 

7.071 

6 

wm 

2 

0.005 

0.0005 

5.0 

7.071 

11 

12 

2 

0.005 

0.00025 

5.0 

7.071 

21 

Table  7:  FMP2DT  Run  Summary  for  TWIGL  Comparison 


For  run  2,  FMP2DT  divided  this  problem  into  two  intervals.  The  first  went  from 
t  =:0tot  =  0.005  sec,  and  the  second  went  from  t  =  0.005  to  0.01  sec.  However,  for 
the  first  interval  only  one  time  step,  with  a.  At  =  0.005  sec,  was  used.  The  second 
interval  used  five  time  steps  with  At  =  0.001  sec.  Thus,  run  2  used  a  total  of  six 
time  steps.  The  answers  for  runs  1  and  2  were  exactly  the  same. 
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This  demonstrates  the  ability  for  FMP2DT  to  yield  a  substantial  savings  in 
computational  cost  due  to  the  fact  that  it  has  an  implicit,  numerical  stable,  time 
discretization.  FMP2DT  demonstrated  this  same  characteristic  for  the  other  com¬ 
binations  of  Ax,  Az  and  At  configurations  shown  in  Table  7. 

The  dilferences  in  the  peak  magnitudes  for  the  TWIGL  and  FMP2DT  calcu¬ 
lations,  shown  in  Figure  14,  yield  no  special  concern  since  the  large  mesh  spacing 
used  makes  the  accuracy  of  either  calculation  questionable.  However,  since  the  flux 
shapes  are  similar,  it  can  be  said  that  the  source  vectors  for  the  two  codes  were 
similar.  This  is  significant  because  both  codes  modeled  a  precursor  source. 

We  benchmark  FMP2DT  now  using  exact  flux  values.  This  will  not  only  vali¬ 
date  the  ability  of  FMP2DT  to  obtain  the  correct  flux  shape,  but  will  enhance  its 
credibility  foi  computational  accuracy. 
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7.2  Slab  Geometry  Benchmark 

Since  the  TWIGL  mesh  spacing  was  much  larger  than  a  MFP,  the  most  that  can 
be  said  about  that  comparison  is  that  FMP2DT  produced  a  similarly  shaped  flux. 
The  TWIGL  work  then  needs  to  be  supplemented  with  some  exact  calculations  to 
establish  FMP2DT’s  ability  to  produce  credible  results. 

B.  D.  Ganapol  [Ref.  29]  published  a  paper  giving  exact  results  for  infinite  one 
dimensional  slab  and  cylindrical  geometries.  FMP2DT  modeled  these  configurations 
by  making  the  mediums  so  large  that  neutrons  born  because  of  a  pulsed  source  at 
time  i  =  0  sec  did  not  have  time  to  leak  out  of  the  medium.  Ganapol  tabulated 
these  calculations  to  aid  in  debugging  programming  errors.  This  section  compares 
FMP2DT  calculations  with  GanapoPs  exact  calculations  for  XZ  slab  geometry. 

Table  8  shows  the  exact  flux  in  an  infinite  slab  with  an  isotropic  pulsed  plane 
source  at  x  =  0  in  a  nonabsorbing  medium.  For  each  time  step,  the  flux  is  calculated 
at  the  mean  free  paths  shown.  The  solutions  were  generated  using  Neumann  series 
for  the  angular  and  scalar  fluxes.  The  neutron  velocity,  u,  was  set  to  be  1  cm/sec, 
and  the  total  cross  section,  S^,  was  set  at  1.0  cm"^.  (This  was  the  case  for  both 
XZ  and  RZ  geometry.)  The  media  for  both  of  the  infinite  geometries  were  non¬ 
multiplying.  Table  9  shows  the  FMP2DT  results,  and  Table  10  shows  the  percent 
relative  error. 

The  data  in  Table  9  were  calculated  w’ith  a  two  dimensional  slab  configuration 
with  a  reflective  boundary  at  z  =  0  cm,  vacuum  boundaries  at  x  =  45  cm  and 
2  =  0  and  2  =  46  cm.  These  values  represent  the  physical  dimensions  of  the  slab 
modeled.  The  fluxes  were  calculated  at  x  =  1, 2,  3,  4,  5,  and  6  cm  respectively  with 
a  corresponding  value  of  2  =  23.0  cm.  These  values  of  x  were  1  MFP  apart,  with 
the  first  being  1  MFP  from  the  boundary.  For  all  runs  in  both  slab  and  cylindrical 


geometries,  the  source  had  a  thickness  of  6.25E-02  cm-  To  insure  convergence,  ten 
runs  were  done  with  different  mesh  spacing  and  different  time  steps.  As  U:e  mesh 
spacing  became  finer,  it  became  necessary  Co  model  the  infinite  medium  by  using 
reflective  boundary  conditions  on.  the  top  and  bottom  of  the  slab,  which  left  only 
the  right  side  allowing  any  leakage.  It  was  found  that  only  five  mesh  points  in  the 
z  direction  were  necessary.  Most  runs  :yere  modeled  with  a  reflector  at  z  =  0  and 
z  =  4  cm  and  the  fliix  data  calculated  at  z  =  2  cm.  In  the  Ccises  where  Az=2  cm, 
the  reflectors  were  put  at  z  =  0  and  z  =  8  cm  with  the  flux  calculations  made 
at  z  =  4  cm.  The  reflective  top  and  bottom  boundary  configurations  yielded  the 
same  answers  as  did  the  configuration  with  vacuum  boundary  conditions  on  the 
top  and  bottom.  Since  the  reflective  configurations  had  fewer  axial  mesh  points, 
there  was  a  substantial  savings  in  computational  cost,  and  they  had  much  fcister.run 
times.  These  same  schemes  were  done  with  the  cylindrical  calculations.  Table  11 
summarizes  the  different  configurations. 

The  results  of  the  FMP2DT  flux  are  extremely  good.  The  error  at  1  mean 
free  path  (MFP)  is  most  noticeable  at  early  time.  It  is  expected  that  the  flux  at 
1  MFP  would  yield  the  most  error  since  this  is  closest  to  the  boundary  and  a  Pi 
approximation  is  more  likely  to  be  suspect  there.  Also,  while  mathematically  the 
line  source  can  be  turned  on  and  off  at  a  time  t  =  0,  at  a:  =  0  cm,  using  a  delta 
function,  the  code  cannot.  The  source  had  to  have  some  finite  dimension,  and  it 
had  to  be  left  on  at  some  finite  time.  These  dimensions  were  extremely  small,  but 
produced  some  error. 

Figure  16  shows  that  the  error  dies  out  as  time  increases.  However,  even  for 
earlier  times,  the  flux  shape  is  in  good  agreement.  Figure  17  shows  the  comparison 
between  the  FMP2DT  and  the  exact  flux  two  MFPs  from  the  source.  Again,  the 
early  time  values  have  the  worst  error.  Figures  18, 19, 20,  and  21  show  the  FMP2DT 
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and  exact  flux  calculations  at  3,  4,  5,  and  6  MFPs  respectively.  The  error  for  these 
calculations  is  extremely  small,  especially  after  the  first  few  time  points.  This  h 
consistent  with  the  expected  flux  behavior  using  a  Pi  approximation. 


TIME 

1  MFP 

2  MFPs 

3  MFPs 

4  MFPs 

5  MFPs 

6  MFPs 

0 

O.OOOOE+00 

OjOOOOE4-00 

O.OOOOE-FCO 

O.OQO0E4-0O 

0.0000E44X) 

0.000064-00 

1 

1.839-lE-Ol  . 

O.OOOOE4-00 

O.OOOOS+00 

O.OOOOE-FOO 

OiXX)OE-FOO 

0.000064-00 

3 

23342E-01 

93836E-02 

8.2978E-03 

O.OOOOE+00 

0.0000E4-00 

0.000064-00 

5 

1.9957B-01 

1.2105E-01 

4.9595E-02 

1.1823E-02 

6.7379E4M 

0.000064-00 

T 

1.7347E-0I 

1.2293E-01 

6.802SE-02 

2a447E-02 

8.415SE4)3 

1303664)3 

9 

1.552SB4)1 

1.1935E-01 

4J01S6E-02 

1.7004E^ 

537«64)3 

11 

l.‘1175E-01 

1.1454  E-01 

4.7953E-02 

2.4433&02 

1.011964)2 

13 

iai20E-Ol 

in9S3E-01 

8.1^E-02 

.53024E-02 

3.0372E.02 

1313764)2 

15 

1.2269E-01 

1.0514E-0I 

8.1158E-02 

5.S305E>02 

1337664)2 

17 

1.1564E-01 

1.009SE-01 

8.043SE-02 

5a390E-02 

3.853IE4J2 

2304164)2 

19 

1.096SE-O1 

9.7166E-02 

7.9349E-02 

5i)663B4)2 

4.1241E4)2 

23150602 

21 

l.ajS5B-01 

93719E-02 

7.S0S6E-02 

SJ0377E-02 

4330564)2 

2376164)2 

23 

l.(K)07E-01 

9XI583E-02 

7.6693B-02 

6J0698E-02 

4.4868B-02 

3.094264)2 

25 

8.7720B-02 

6J0744B-02 

4.6042B4)2 

3.2757602 

27 

8.5093E-02 

7aS85&02 

4.691264)2 

3.426564)2 

29 

8iM60E-02 

S.26SSB4» 

6J0301E-02 

4.754?E4)2 

33515602 

31 

S.6606E-02 

Sa4&fE-02 

S:9910E-02 

4.793464)2 

3.654964)2 

33 

8.4009E-02 

6OS72E-02 

53448E-02 

4327464)2 

3.740064)2 

35 

S.I632E-02 

7iM91E-02 

6.8624B-02 

S.S937E-02 

4344564)2 

3309964)2 

37 

6.7424E-02 

53393EM)2 

4351964)2 

3366964)2 

39 

MOaiekesaM 

73041E432 

4351564)2 

33129602 

•il 

7.1479E.02 

6.5I67E-02 

5.7247B-02 

4345064)2 

3349764)2 

•13 

6.410SE-02 

337^602 

45 

mmmpifjm 

6S630E-02 

6309I&02 

5£074E-&2 

4317764)2 

4.000764)2 

Table  8:  Exact  FIu.x  Due  to  an  Isotropic  Pulsed  Plane  Source  at  x  =  0  in  a  Nonab¬ 
sorbing  Infinite  Medium  For  Slab  Geometry 
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TIME 

_  1  MFP 

2MFPs 

3  MFPs 

4  MFPs 

5  MFPs 

6  MFPs 

0 

O.OOOOE+00 

O.OOOOB+00 

O.OOOOE+00 

O.OOOOE+00 

O.OOOOE+00 

O.OOOOE+00 

1 

O.OOOOE+00 

O.OOOOE+00 

O.OOOOB+00 

O.OOOOE+00 

O.OOOOE+00 

3 

2.7485E-01 

7.4742E-02 

O.OOOOE+00 

O.OOOOE+00 

O.OOOOE+00 

5 

2.0536E-01 

1.3069E-01 

4.3453B-02 

9.2030E-03 

1.4897E-03 

O.OOOOE+00 

7 

1.6964E-01 

1.3084E-01 

7.0663B-02 

2.5537E-02 

6.6768E-03 

1.3782E-03 

9 

1.5137B-01 

1.2190E-01 

8.0677E-02 

4.0570E-02. 

1.5307E-02 

4.4989E-03 

11 

1.3853B-01 

1.1488E-01 

8.2818E-02 

4.9817E-02 

2.4089E-02 

9.3242E-03 

13 

1.2850E-01 

1.0927E-01 

8.2747E-02 

5.4848E-02 

3.0971E-02 

1.4612E-02 

15 

1.2038E-01 

1.0442&01 

8.1953B-02 

5.7699E-02 

3.5895E-02 

1.9389E-02 

17 

1.1362E-01 

1.0012E-01 

8.0792E-02 

5.9390E-02 

3.9427E-02 

2.3361E-02 

19 

1.0789B-01 

9.6275E-02 

7.9427B-02 

6.0360E-02 

4.2021EJ-02 

2.6594E-02 

21 

1.0294E-01 

9.2818E-02 

7.7963B-02 

6.0845E-02 

4.3951E-02 

2.9232E-02 

23 

9.8615B-02 

8.9692E-02 

_  7.6467E-02 

6.0993B-02 

4.5387B.02 

3.1392E-02 

25 

9.4792E-02 

8.6850E-02 

7.4978B.02 

6.0904E-02 

4.6448E-02 

3.3166E>-02 

27 

9.1382E-02 

8.42S4E-02 

7.3519E-02 

6.0649B-02 

4.7220E-02 

3.4624E-02 

29 

8.831SE-02 

8.1871B-02 

7.2104B-02 

4.7766&02 

3.5822B-02 

31 

8.5538B-02 

7.9674B-02 

7.0738E-02 

5.9821E-02 

4.8136B.02 

3.6805E-02 

33 

8.3006E-02 

7.7642E-02 

6.9425B-02 

5.9309E-02 

4.8366I5-02 

3.7609&02 

35 

8.0687E-02 

7.5755E-C2 

6.8166E-02 

5.8758E-02 

4.8485E-02 

3.8263&02 

37 

7.8552E-02 

7.3997B-02 

6.6960B-02 

5.8182E-02 

4.8515E-02 

3.8792B-02 

39 

7.6578B-02 

7.2355E-02 

5.7590E-02 

4.8475E-02 

3.9216E-02 

41 

7.4745B-02 

7.0815E-02 

5.6991E-02 

4.8377E-02 

3.9552E-02 

43 

7.3038E-02 

6.9369B-02 

6.3642B-02 

5.6389E-02 

4.8234E-02 

3.9811E-02 

45 

7.1442B-02 

6.8007B-02 

6.2629E-02 

5.5789E-02 

4.8054B-02 

4.0007E-02 

Table  9:  FMP2DT  Flux  Due  to  an  Isotropic  Pulsed  Plane  Source  at  a:  =  0  in  a 
Nonabsorbing  Medium  For  Slab  Geometry 
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TIME 

1  MFP 

2  MFPs 

3  MFPs 

4  MFPs 

5  MFPs 

6  MFPs 

1  ; 

43.13 

N/A 

N/A 

N/A 

N/A 

N/A 

3 

14.80 

20.35 

19.14 

N/A 

N/A 

N/A 

5 

2.90 

7.96 

12.38 

22.16  : 

121.09 

N/A 

7 

2.21 

6.43 

3.87 

10.23 

20.66 

8.34 

9 

2.52 

2;i4 

5.62 

0.96 

9.98 

19.32 

11 

2.27 

0.30- 

3.54 

3.89 

1.41 

10.51 

13 

2.06 

0.38 

1.91 

■sai 

1.97 

3.47 

15 

1.88 

0.68 

0.98 

2.94 

2.60 

0.07 

17 

1.75 

0.83 

0.44 

1.71 

2.33 

1.39 

19 

1.63 

0.92 

0.10 

1.22 

1.89 

1.70 

mm 

1.54 

0.96 

0.13 

0.78  : 

1.49 

1.64 

23 

1.45 

0.98 

0.29 

0.49 

1.16 

1.45 

25 

1.39 

0.99 

0.41 

0.26 

0.88 

1.25 

27 

1.33 

0.99 

0.50 

0.09 

0.66 

1.05 

29 

0.99 

0.56 

0.04 

0.47 

0.86 

31 

1.23 

0.98 

0.60 

0.15 

0.32 

0.70 

33 

1.19 

0.97 

0.64 

0.23 

0.19 

0.56 

35 

1.16 

0.96 

0.67 

0.30 

0.08 

0.43 

37 

1.13 

0.95 

0.69 

0.36 

0.01 

0.32 

39 

1.10 

0.94 

0.70 

0.41 

0.08 

0.22 

41 

1.07 

0.93 

0.72 

0.45 

0.15 

0.14 

43 

1.05 

0.92 

0.73 

0.48 

0.21 

0.06 

45 

1.03 

0.91 

0.73 

0.51 

0.26 

0.00 

Table  10:  Tlie  Percent  Relative  Error  for  the  FMP2DT  Slab  Calculations 


run 

reflectors 

Ax  or  Ar,  Az 

At 

1 

no 

1.0  cm 

1.0  sec 

2 

yes 

0.5  cm 

1.0  sec 

3 

yes 

0.5  cm 

0.5  sec 

4 

yes 

0.5  cm 

0.25  sec 

5 

yes 

1.0  cm 

2.0  sec 

6 

yes 

1.0  cm 

1.0  sec 

B 

yes 

1.0  cm 

0.5  sec 

8 

yes 

2.0  cm 

4.0  sec 

9 

yes 

2.0  cm 

2.0  sec 

10 

yes 

2.0  cm 

1.0  sec 

Table  11:  Run  Summary  for  the  Infinite  Slab  and  Cylindrical  Geometries 
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(08S/3**ujo/u)  xny  uoJinaN 

Figure  16:  FMP2DT  and  Exact  Flux,  XZ  Geoinetry  1  MFP  From  Boundary 
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8.0  2 

Time  (sec 


{O0s/3*»ujo/u)  xny  uoJinaN 

Figure  17:  FMP2DT  and  Exact  Flux,  XZ  Geometry  2  MFPs  From  Boundary 
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Figure  18:  FMP2DT  and  Exact  Flux,  XZ  Geometry  3  MFPs  From  Boundary 
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Figure  19:  FMP2DT  and  Exact  Flux,  XZ  Geoinetry  4  MFPs  From  Boundary 
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Figure  20:  FMP2DT  and  Exact  Flux,  XZ  Geometry  5  MFPs  From  Boundary 
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0.00065 


{08S/3*#ujo/u)  xny  uoJinaN 

Figure  21:  FMP2DT  and  Exact  Flux,  XZ  Georitetry  6  MFPs  From  Boundary 
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7.3  Infinite  Cylindrical  Geometry  Benchmark 

The  Ganapol  paper  already  cited  also  gave  results  for  cylindrical  geometry.  Ta¬ 
ble  12  shows  the  exact  flux  due  to  an  isotropic  pulsed  line  source  a  r  =  0  cm  in  a 
nonabsorbing  infinite  medium. 

As  with  XZ  geometry,  several  different  runs  were  made  to  insure  convergence. 
The  run  summary  is  given  in  Table  11.  For  RZ  geometry,  there  is  a  reflective 
boundary  at  the  radial  center,  i.e.  at  r  =  0  cm.  The  rest  of  the  boundary  conditions 
for  the  different  configurations  are  cis  described  for  the  XZ  calculations. 

Again,  just  like  the  XZ  case,  FMP2DT  calculated  the  flux  with  good  agreement 
with  the  exact  flux.  Table  12  shows  the  exact  flux,  Table  13  shows  the  FMP2DT 
calculations,  and  Table  14  shows  the  percent  relative  error  with  respect  to  the  exact 
flux.  As  with  the  XZ  geometry,  the  line  source  had  a  small  finite  thickness,  and 
the  source  was  turned  on  and  off  with  some  small  finite  time.  This  introduces  some 
error  automatically  that  cannot  be  omitted. 

Figures  22,  23,  24,  2.5,  26,  and  27  show  the  graphs  comparing  the  exact  and 
FMP2DT  fluxes  at  1  MFP,  2  MFPs,  3  MFPs,  4  MFPs,  5  MFPs,  and  6  MFPs 
respectively.  As  in  the  XZ  geometry  calculations,  the  worst  error  occurred  at  the 
earlier  times.  This  is  caused  by  the  characteristic  of  the  Pi  approximation  and  the 
finite  source  configuration  already  stated.  .As  the  calculations  moved  away  from 
the  cylinder’s  center  line,  the  F.v'IP2DT  flux  was  nearly  identical  to  the  exact  flux 
reported  by  Ganapol. 
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TIME 

1  MFP 

3  MFPs 

4MFPs 

5  MFPs 

6  MFPs 

1 

O.OOOOE+00 

OJ0000E4-0O 

•O.OOOOEfOO 

O.OOOOEtOO 

O.(XX)0E4-CO 

O.OOOOE4-00 

3 

3.1734 E-02 

O.OOOOE+00 

O.OOOOEfOO 

O.OOOOE+00 

O.OOOOEfOO 

5 

4.5377E-02 

2.8563E-02 

i.22S5B-02 

3.2456E-03 

O.OOOOE4-00 

O.OOOOE+00 

7 

3^54E-02 

23830E-02 

1J458E-02 

58172E-03 

1.8173E-03 

9 

2.6027E-02 

2.0134E-02 

1.3033E4)2 

69753E-C3 

3.0260E-03 

11 

2.1376E-02 

1.7346E-02 

1.2201E4)2 

73928E-03 

.38227E-03 

13 

1^133E-02 

I.5205E-02 

1.1312E-02 

430-l2E4)3 

2.1729E-03 

15 

1.5744E-02 

1.3521 E-02 

1.W76E-02 

7J063E-03 

4.5721 E-03 

2.5552E-03 

17 

l^llB-02 

1.2165E-02 

9.7197E-03 

7.0S42E-03 

4.699SE-03 

19 

1.2460E-02 

1.1053E-02 

9.0465B.03 

68239E-03 

4.7367E-03 

21 

1.1282E-02 

1J)125E.02 

8.4492E-03 

6.5513E-03 

4.7150E4)3 

3.1447E-03 

23 

1.030SB-02 

9J393E-03 

7.9191E-03 

6.2809E-03 

4.65S8E-03 

3.2218E-03 

25 

9.4S93E-03 

S.6659E-03 

7.1470E-03 

6D192E-03 

4.5731E-03 

27 

S.7907B-03 

SD825E-03 

5.7700E-03 

4.4761E-03 

29 

3.187SE-03 

6.6462E-03 

S.5346E-03 

4.3711E4)3 

31 

7.6623E-03 

7.1224E-03 

6.3047E-03 

53135E-03 

4.2620E4» 

33 

6.7228E-03 

5.9955E-03 

5.1062E-03 

4.1S1SE03 

3.221SE-03  ! 

35 

6.7907E-03 

6.3654E-03 

5.7144E-03 

49122E-03 

37 

6.4252E-03 

61M41E-03 

5.4579E-03 

4.7305E-03 

3.9347E-03 

39 

6.0971E.03 

5.7535E-03 

5.2229E.03 

4.S60SE-03 

3.S29SE4)3 

41 

5.800SB-03 

5.4S95E-03 

5.0071B-03 

3.72S1E-03 

43 

5^20E-03 

5.24S7E-03 

4.S0S0E-03 

4.2520E-03 

3.6299E-03 

45 

5.2870E-03 

5D2S0E-03 

4.6240E-03 

4.1I19E-03 

3.5352E4)3 

Table  12;  Exact  Flux  Due  to  an  Isotropic  Pulsed  Plane  Source  at  r  =  0  in  a 
Nonabsorbing  Infinite  Medium  For  RZ  Geometrj' 


TIME 

■DSQIH 

2  MFPs 

3  MFPs 

4  MFPs 

5  MFPs 

6  MFPs 

1 

O.OOOOEfOO 

O.OOOOE+00 

O.OOOOE+00 

O.OOOOE4-CO 

O.OOOOE-100 

3 

9.7929E-02 

2.7334 E-02 

O.O00OE4-OO 

OXOCOE4-00 

O.OOOOE4-00 

O.OOOOEfOO 

5 

33S92E-02 

1.1493E-02 

2.3169E-03 

O.OOOOE+00 

O.0000E4-00 

7 

2.6S79E-02 

2.5258E-02 

1.4980E-02 

5.5095E-03 

1.4101E-03 

9 

2.3720E-02 

li)766E-02 

I.4075E-02 

7.4456E-03 

2.S511E-03 

11 

1.9968E-02 

1.6787E-02 

1.2558E-02 

3.9532E4)3 

1.54g2E-03  1 

13 

1.7160E-02 

1.-172VE-02 

I-1371E-02 

4.5339E-03 

15 

1.S027E-02 

1.3124  E-02 

1.0431E-02 

4.7768E-03 

2.6427B-03  1 

17 

1.3360E-02 

1.1832E-02 

9.6363B-03 

7.1S07E-03 

4.8536E4)3 

19 

1.2023E-02 

1D770E-02 

S.9iS5E-03 

68702E-03 

4AI72E-03 

3.1in)E-03  1 

21 

1.0D28E-02 

9882SE-03 

8.3475E-03 

6.5G71E-03 

4.7933E-03 

23 

1.0015E-02 

9.1301  B-03 

7.S1SSE-03 

6.2770E-03 

4.7I02E4)3 

25 

9.2!I9E-03 

8.4S36E-03 

BSiSmm 

6D029E-03 

4.6097E-03 

27 

S.579-IE-03 

6.9312E-03 

5.7457E-03 

4.499-iE-03 

29 

8.0052E-03 

6.56I0E-03 

5.5055E-03 

4JS-I3E-03 

3.3077E-03 

31 

65961E-03 

6.2251E-03 

5.2815E-03 

4.2677E-03 

3.2799E-03 

33 

6i»95E-03 

5.92I3E-03 

5D727B-03 

3.2429E-03 

35 

6.6C61E-03 

e.263iE-03 

S.&i52E-03 

48760B-03 

4.0331  E4)3 

3.I993E-C3 

37 

G.3I38E-03 

5.95I6E-03 

5J93-1E-03 

3.9272E-03 

3.1511E-03 

39 

5.9969E-03 

5.6694  E-03 

S.1628E-03 

4.52«)E-03 

3.8I99B4)3 

41 

S.7I02B-03 

5.<:i26E-03 

43681E-03 

3.716-IE-03 

3.OI6SE.03 

43 

5.4I06E-03 

5.17SOE-03 

■m-fra.y 

4.219DE-03 

3.G169E-03 

2.9927E-03 

45 

5.211SE-03 

4.9629E-03 

4.5747E-03 

3.5213E4)3 

Table  13:  FMP2DT  Flux  Due  to  an  Isotropic  Pulsed  Plane  Source  at  r  =  0  in  a 
Nonabsorbing  Infinite  Medium  For  RZ  Geometry 


S5 


TIME 

1  MFP 

2  MFPs 

3  MFPs 

4  MFPs 

5  MFPs 

6  MFPs 

3 

32.36 

13.87 

N/A 

N/A 

N/A 

N/A 

5 

9.81 

18.66 

6.43 

28.61 

N/A 

N/A 

7 

13.16 

5.99 

11.31 

5.29 

22.41 

22.14 

9 

8.86 

1.83 

7.76 

6.74 

5.78 

19.31 

322 

2.93 

7.05 

3.41 

6.83 

3.14 

0.52 

4.47 

5.34 

0.73 

15 

4.55 

2.94 

0.43 

2.51 

4.48 

3.42 

17 

3.96 

2.74 

0.86 

1.36 

3.27 

3.71 

19 

3.51 

2.56 

1.08 

0.68 

2.33 

3.24 

21 

3.14 

2.39 

1.20 

0.24 

1.66 

2.66 

23 

2.84 

2.24 

1.27 

0.06 

1.17 

2.15 

25 

2.61 

2.10 

1.29 

0.27 

0.80 

1.72 

27 

2.40 

1.98 

1.29 

0.42 

0.52 

1.37 

29 

2.23 

■mi 

1.-28 

0.53 

1.08 

31 

2.08 

1.77 

1.26 

0.60 

0.13 

O.Sd 

33 

1.95 

1.69 

1.24 

0.66 

~0.0 

0.65 

35 

1.83 

1.60 

1.21 

0.70 

0.10 

0.49 

37 

1.73 

1.53 

1.18 

0.72 

0.19 

0.36 

0.36 

0.07 

0.39 

0.00 

Table  14:  The  Percent  Relative  Error  for  the  FMP2DT  Cylindrical  Calculations 
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Figure  22:  FMP2DT  and  Exact  Flux,  RZ  Geometry  1  MFP  From  Boundary 
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Figure  23:  FMP2DT  and  Exact  Flux,  RZ  Geometry  2  MFPs  From  Boundary 


88 


(oas/3#*uJO/u)  xny  uoJinaN 

Figure  24;  FMP2DT  and  Exact  Flux,  RZ  Geometry  3  MFPs  From  Boundary 
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Figure  25:  FMP2DT  and  Exact  Flux,  RZ  Geometry  4  MFPs  From  Boundary 
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Figure  26:  FMP2DT  and  Exact  Flux,  RZ  Geornetry  5  MFPs  From  Boundary 
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Figure  27;  FMP2DT  and  Exact  Flux,  RZ  Geometry  6  MFPs  From  Boundary 
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7.4  Benchmark  Conclusions 


FMP2DT  could  not  be  benchmarked  with  TWIGL.  The  TWIGL  calculations  were 
done  with  a  much  too  coarse  a  mesh  spacing.  However,  comparing  FMP2DT  with 
TWIGL’s  flux  shape  showed  that  the  precursor  contribution  was  at  least  generating 
a  similar  source  vector  as  was  TWIGL’s  source  vector.  And  by  presenting  the 
TWIGL  comparison,  a  demonstration  of  some  of  FMP2DT’s  input  characteristics 
was  accomplished.  The  benchmarks  with  the  exact  flux  in  infinite  XZ  and  RZ 
geometries  did  establish  FMP2DT’s  computation  creditability. 

Since  a  Pi  approximation  is  most  suspect  near  boundaries  and  strong  sources,  it 
is  most  likely  that  FMP2DT  would  show  the  greatest  error  nearest  to  the  boundary 
or  source.  Figures  16  and  17  show  this  to  be  true.  Angular  flux  approximations 
usually  do  not  predict  the  early  time  behavior  of  the  flux  because  the  wave  behavior 
of  the  flux  is  not  fully  accounted  for  [Ref.  29].  However,  given  time  after  the  source 
is  turned  off,  the  error  disappears,  even  in  close  proximity  to  the  boundary.  This 
implies  that  the  flux  has  time  to  become  more  isotropic  so  that  the  Pi  approximation 
is  a  better  representation  of  the  angular  flux. 

The  fact  that  a  Pi  approximation  has  difficulty  predicting  the  early  behavior 
of  the  flux  is  compounded  by  the  fact  that  the  source  was  0.0625  cm  thick.  This 
means  that  at  one  second,  the  leading  edge  of  the  wave  is  at  1.0625  cm  because  the 
neutron  velocity  is  1  cm/sec.  Likewise,  there  is  an  offset  in  the  leading  edge  at  3  and 
5  sec.  This  implies  that,  at  these  times,  the  corresponding  calculations  at  1,  3,  and 
5  MFPs  are  not  predicting  the  leading  edge  of  the  wave  since  it  has  already  passed. 
Looking  at  the  percent  relative  error  for  slab  geometry  in  Table  10,  the  largest  error 
for  these  MFPs  is  at  times  <  =  1,  i  =  3,  and  i  =  5  sec.  Looking  at  Tables  8  and  9, 
at  3  and  5  MFPs,  it  appears  that  the  flux  is  over  predicted.  This  supports  the  fact 
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that  the  leading  edge  of  the  wave  has  already  pass  through. 

For  RZ  geometry,  similar  behavior  is  true.  There  is  an  inherent  radial  dependency 
for  cylindrical  geometry  that  is  more  pronounced  than  for  slab  geometry.  The  largest 
errors  are  for  the  earlier  time  fluxes.  Looking  at  Figures  23  and  24,  the  early  time 
FMP2DT  flux  shape  at  2  MFPs  appear  to  be  taking  the  shape  of  the  3  MFP  flux 
shape.  This  again  supports  the  suggestion  that  the  leading  edge  of  the  wave  has 
already  passed  through.  Experiments  with  the  mesh  spacing  also  seemed  to  confirm 
this. 

Overall,  the  benchmark  results  with  the  exact  flux  are  remarkably  good.  Usually, 
comparisons  between  transport  codes  and  exact  answers  given  in  literature  are  for 
several  MFPs  from  the  source  and  boundary.  This  is  because  most  codes  are  diffusion 
codes,  and  diffusion  theory  breaks  down  near  the  boundary.  Diffusion  codes  do  not 
represent  the  anisotropic  nature  of  the  leading  edge  of  the  wave  at  early  time  steps. 
After  a  few  MFPs,  the  flux  becomes  more  isotrcoic  because  of  resulting  collisions, 
and  diffusion  theory  becomes  more  valid.  Normally,  a  Pi  approximation  gives  a 
better  representation  of  the  anisotropic  nature  of  the  flux  and  is  better  representation 
of  the  leading  edge  of  the  wave  at  earlier  time.  If  more  accuracy  is  desired  at  early 
time,  then  a  larger  P„  approximation  would  be  a  better  model  for  the  angular  flux. 
But  as  the  jP„  approximations  are  made  and  programmed,  the  computational  costs 
rise  substantially.  Thus,  the  Pi  approximaUv^n  in  FMP2DT  still  is  economically 
attractive,  and  the  early  time  errors  may  be  accept^,^»c  for  most  computations. 

In  the  next  section,  we  examine  some  appljt  .tvions  for  FMP2DT,  and  identify 
some  few  group  problems  where  FMP2DT  can  be  used  as  an  analytical  tool. 
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8  FMP2DT  Applications 


FMP2DT  has  applications  to  many  nuclear  engineering  problems  of  interest  to  a 
number  of  communities.  For  examples,  two  problems  were  chosen,  both  in  RZ  ge¬ 
ometry.  The  first  one  can  only  be  resolvfj.3  decay  of  a  critical  reactor  after 

injection  of  a  pulse  of  neutrons.  This  y,  .-n..  .  Iten  used  to  determine  neutron 
lifetime  estimates  in  a  new  or  research  react*  ..  The  second  experiment  is  represen¬ 
tative  of  most  pulsed  logging  problems.  Tw  '  ..iiensional  effects  can  dominate  the 
measurements  to  the  extent  that  simpler  on^  .-.-nensional  models  are  useless.  These 
problems  are  only  to  demonstrate  FMP2DT.  Therefore,  no  detailed  physics  will  be 
developed.  For  example,  for  the  first  sample  problem,  FMP2DT  will  show  how  it 
can  aid  in  collecting  flux  data  for  a  Rossi-alpha  experiment.  However,  the  actual 
value  of  alpha  will  not  be  calculated. 

8.1  AGN-201  Rossi-Alpha  Problem 

The  first  example  problem  chosen  entailed  modeling  the  UNM  AGN-201  reactor. 
The  AGN-201  is  a  low-power  thermal  reactor  used  mostly  for  training  purposes.  It 
is  a  right  cylinder  core,  25  cm  in  diameter.  The  fuel  is  a  mixture  of  20%  enriched 
UO2  and  polyethylene.  The  critical  mciss  of  the  reactor  is  about  665  grams  of 
The  experiment  modeled  involved  a  pulsed  source  located  ai  the  center  of  the 
UNM  AGN-201  reactor  that  was  initially  operating  at  a  critical  state.  This  type  of 
experiment  is  known  as  a  Rossi-alpha  experiment.  A  pulse  of  neutrons  is  injected 
into  the  reactor  with  a  spatial  distribution  that  excites  multiple  spatial  modes.  After 
a  period  of  time,  all  modes  other  than  the  fundamental  mode  die  out.  The  decay 
of  the  fundamental  mode  then  can  be  related  to  the  neutron  lifetinie  01  generation 
time  in  the  critical  reactor.  The  crucial  piece  of  information  that  the  calculation  can 


provide  is  the  time  at  which  counters  can  be  gated  on  to  measure  the  fundamental 
mode  decay.  If  the  counters  are  gated  on  too  soon,  the  fundamental  mode  will  be 
contaminated  by  measurements  of  the  higher  modes.  If  the  counters  are  gated  on 
too  late,  the  data  will  suffer  from  poor  resolution. 

The  two  dimensional  time  dependent  calculation  can  be  recorded  at  multiple 
locations.  When  the  time  history  of  the  flux  at  all  spatial  locations  has  the  same 
exponen^'  d  decay  on  a  semi-log  plot,  all  higher  modes  have  decayed  out.  At  this 
time  the  counters  can  be  gated  on  and  the  decay  data  recorded. 

8.1.1  AGN  Calculation  Description 

The  geometry  for  this  problem  is  shown  in  Figure  28.  The  source  was  turned  on 
for  one  microsecond  and  then  turned  off.  The  source  energy  was  set  at  0.01  MeV. 
The  magnitude  of  the  oource  energy  was  chosen  arbitrarily  and  poses  no  significance 
here.  The  mesh  was  about  1  cm  in  the  core  for  the  radial  direction,  and  1  cm  in 
the  axial  direction.  The  At  for  th’  problem  was  20  microseconds  for  50  steps.  The 
calculations  were  stopped  at  1.001  millisecond.  The  initial  flux  was  generated  b>  a 
FEMP2D  calculation.  The  data  were  taken  at  the  axial  center  line  which  is  at  z=0 
cm  in  Figure  28.  The  radial  values  were  at  r=0.985,  2.954,  4.923,  6.982,  and  8.861 
cm.  Tliese  radial  positions  are  all  in  the  core.  The  last  position  was  chosen  to  stay 
away  fr*.xn  the  core  and  graphite  interface  since  the  reflection  of  the  neutrons  from 
the  graphite  changes  the  slope  of  the  flux  shape.  A  three  group  problem  was  modeled 
using  cross  section  data  from  an  AMPX  library.  This  was  because  of  the  availability 
of  the  cross  section  set,  and  three  groups  are  enough  to  demonstrate  FMP2DT’s 
multigroup  computational  ability.  Table  15  shows  the  energy  group  structure  for 
this  problem. 
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Group 

E  Max  (eV) 

E  Min  (eV) 

1 

1.7330E-b07 

l.OOOOE-l-07 

2 

l.OOOOE-1-07 

9.9990E-02 

3 

9.9990E-02 

l.OOOOE-05 

Table  15:  Three  Group  Energy  Structure 
8.1.2  AGN  Results 

Figure  29  shows  the  plot  of  the  thermal  flux  at  the  stated  core  positions.  From  the 
graph,  it  can  be  seen  that  after  about  0.27  milliseconds,  the  slopes  of  the  thermal  flux 
are  the  same  for  all  the  radial  positions.  This  implies  that  all  but  the  fundamental 
flux  mode  has  died  out.  This  is  the  time  to  gate  the  counters  on  and  record  the  flux 
data.  Any  time  previous  would  result  in  higher  mode  data  contamination  since  the 
slope  is  not  straight  for  all  spatial  points. 

Next,  we  turn  our  attention  to  a  uranium  logging  problem.  It  also  entails  ob¬ 
serving  neutron  flux  decay  much  like  the  Rossi-alpha  observations  done  here. 
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Figure  29:  Flux  Decay  for  the  AGN-201  Reactor 
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8.2  Uranium  Logging  Problem 

The  search  for  uranium  in  soil  is  accomplished  by  what  is  called  a  uranium  borehole 
logging  procedure.  A  large  hole  is  bored  into  ore-bearing  rock  formations,  and  either 
a  pulsable  source  or  some  neutron  source  such  as  ^®^Cf  is  inserted.  The  resulting 
decay  of  prompt  fission  neutrons  in  epithermal  energy  region  is  observed.  Rock 
formations  with  uranium  present  will  show  an  increase  for  short  time  as  a  result 
of  fission  events,  wherccis  formations  without  uranium  will  show  a  constant  rate  of 
epithermal  neutron  population  decay. 

The  model  chosen  was  the  prompt  neutron  logging  problem  developed  by  James 
H.  Renken  (Ref.  32].  Renken’s  work  was  performed  in  one  dimensional  geometry 
and  assumed  that  a  14  MeV  source  and  a  neutron  detector  were  co-bcated.  This  is 
physically  impossible.  In  fact  they  are  separated  by  lO’s  of  centimeters.  Since  this  is 
true,  diffusion  of  the  thermal  neutrons  must  be  considered  in  the  analysis  as  well  as 
absorption.  A  two  dimensional  model  provides  a  much  better  tool  to  analyze  data 
from  this  physical  problem. 

When  the  logging  problem  is  modeled,  the  fundamental  problem  that  is  corre¬ 
lated  with  the  material  properties  of  the  interrogated  rock  is  the  long  time  (hundreds 
of  microseconds)  decay  constant  of  the  thermal  neutron  population.  This  decay  con¬ 
stant  will  vary  with  distance  from  the  source  if  thermal  neutron  diffusion  is  present. 
The  analysis  presented  here  demonstrates  this  to  be  the  case  and  indicates  that  an 
individual  probe  must  be  calibrated  for  the  source  to  detector  spacing  involved.  This 
type  of  analysis  can  be  useful  for  all  pulsed  logging  systems. 


100 


8.2.1  Logging  Calculation  Description 


The  logging  geometry  is  shown  in  Figure  .30.  .4  borehole  is  drilled  into  potential  ore 
bearing  rock.  The  borehole  has  a  diameter  of  about  4  cm  and  is  filled  with  water. 
A  probe,  about  3  cm  in  diameter  and  60  cm  in  length,  is  inserted  into  the  borehole. 
The  detector  is  about  4  cm  from  the  source.  The  source  is  about  1  cm  in  diameter 
and  is  located  in  the  center  of  the  probe  as  shown. 

The  same  three  group  stru-ture  as  used  for  the  previous  example  is  used  for  this 
problem.  The  group  structure  is  shown  in  Table  lo.  The  radial  mesh  was  about  1  cm 
in  the  axial  and  radial  directions.  (Except  where  material  interfaces  occurred.)  The 
source  energy  was  14  MeV,  and  it  was  turned  on  for  10  niicroseconds.  The  flux  data 
were  then  tabulated  for  -590  microseconds  with  a  A£  =11.8  microseconds  for  50  time 
steps.  The  soil  contains  a  mixture  of  ^U,  ^U,  H2O,  and  SiOj.  The  production  of 
delayed  neutrons  is  neglected  because  their  numbers  will  be  small  during  the  short 
counting  period.  A  reflective  boundary  condition  is  used  at  r  =  0  cm  and  vacuum 
boundary  conditions  are  used  at  r  =  35  cm,  z  =  0  cm,  and  s  =  60  cm. 

S.2.2  Logging  Problem  Results 

Figure  31  shows  the  graph  of  the  flux  decay  at  r  =  0  cm  and  axial  positions  of 
r  =  38,  44,  50,  56,  and  59  cm.  The  source,  also  located  at  r  =  0  cm,  was  between 
z  =  Zl  and  32  cm.  .4s  shown,  the  flux  has  diflerent  slopes  for  cacli  axial  position. 
This  diflerent  slope  data  clearly  show  that  cacli  position  docs  have  a  diflerent  decay 
constant.  This  must  be  taken  into  account  in  detector  calibration  lo  insure  credible 
counting. 
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9  Conclusions 


The  new  code,  FMP2DT,  demonstrated  most  of  the  characteristics  that  were  desired 
before  this  research  began.  Indeed,  its  benchmarking  turned  out  to  be  better  than 
expected.  Its  input  deck,  while  somewhat  tedious  for  source  transients  with  many 
time  steps,  is  fairly  straight  forward  with  respect  to  its  sister  FEMP  codes.  The 
following  is  to  identify  some  of  FMP2DT’s  characteristics. 

The  nature  of  having  to  calculate  the  radial  and  axial  currents,  ’9?io  and  ^n, 
means  having  to  invert  two  matrices.  Depending  on  the  mesh  size  of  the  problem, 
and  the  number  of  energy  groups  considered,  this  can  be  a  major,  and  very  costly, 
operation.  Also,  the  Calculation  of  the  precursor  concentrations,  C,  can  be  even 
more  computationally  expensive  depending  on  the  number  of  precursor  families  to 
be  modeled.  All  of  these  computations  are  magnified  by  the  number  of  time  steps 
chosen.  Therefore,  selecting  the  spatial  mesh  wisely,  and  using  a  group  structure 
that  is  just  fine  enough  to  satisfy  the  problem  physics  is  prudent.  However,  the  user 
may  not  know  these  optimum  configurations.  In  that  event,  a  guess  for  an  equal 
Aa;  and  may  be  selected  and  a  Ai  chosen  to  give  a  tolerable  computational 
time.  Then  some  variations  from  these  values  can  be  used  to  insure  convergence. 
Determination  of  the  proper  group  structure  is  subject  to  choosing  the  number  of 
groups,  running  the  problem,  collapsing  the  groups,  and  observing  the  flux  change. 

As  stated  in  earlier  sections,  the  Euler  backward  differencing  scheme  for  time  dis¬ 
cretization  is  inherently  stable.  However,  the  user  should  be  careful  not  to  assume 
convergence  upon  an  initial  run;  stability  does  not  imply  convergence.  This  is  espe¬ 
cially  true  for  time  dependent  problems,  FMP2DT’s  answer  needs  to  be  compared 
with  more  than  one  run  for  any  new  problem. 

The  ability  for  a  calculation  of  xf  is  unique.  These  values  are  very  difficult  to 


find.  However,  if  data  are  available  to  establish  the  delayed  group  spectrum,  then 
it  can  be  entered  in  FMP2DT’s  input  deck.  In  fact,  if  any  cross  section  parameters 
are  available,  those  may  be  entered  likewise. 

The  benchmarking  of  EMP2DT  established  FMP2DT’s  computational  accuracy. 
The  errors  shown  in  the  early  times  for  both  XZ  and  RZ  geometries  are  consistent 
with  theory.  With  differing  neutron  energy  groups,  where  the  neutron  velocities  are 
much  faster,  the  early  time  error  may  not  be  as  pronounced. 

The  few  group  problems  demonstrated  some  practical  applications  for  FMP2DT. 
The  AGN-201  problem  demonstrated  how  an  experiment,  such  as  Rossi-alpha,  would 
be  accomplished  by  gating  the  counter  in  the  proper  time  interval  to  insure  credible 
data  recording. 

The  uranium  logging  data  showed  that  for  different  spatial  points  in  the  borehole, 
the  flux  decays  at  different  rates.  Therefore,  for  the  detector  to  yield  valuable 
information,  it  must  be  calibrated  for  the  source  to  detector  spacing.  Thus,  the  two 
dimensional  difl'usion  effects  need  to  be  accounted  for. 

FMP2DT  can  be  enhanced  in  many  ways.  Here  are  a  few  suggestions  for  future 
work  and  research. 

9.1  Future  Research 

The  first  thing  that  could  enhance  FMP2DT  is  to  program  the  ability  to  make  the 
cross  sections  temperature  dependent.  This  could  be  done  by  coupling  the  FMP2DT 
equations  with  some  thermal-hydraulic  equations,  and/or  possibly  with  equations 
describing  heat  generation  due  to  particle  interactions.  Included  with  the  particle 
interactions  are  the  effects  of  ganuna  heating.  Of  course,  because  of  computational 
cost  and  complexities,  there  are  limits  as  to  the  detail  for  all  of  this.  But  a:,  least 
some  work  of  this  nature  could  be  feasible. 
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When  doing  the  benchmarks  with  Ganapol’s  data  for  RZ  and  XZ  geometries, 
the  greatest  error  was  at  the  early  times.  A  P„  expansion  of  the  angular  flux  with 
n  >  1  would  be  more  accurate  in  describing  the  angular  flux’s  anisotropic  nature 
at  those  times.  Also,  it  would  permit  a  more  accurate  representation  at  boundaries 
and  material  interfaces.  This  might  be  the  best  enhancement  for  practical  purposes. 

Since  we  have  a  two  dimensional  time  dependent  code,  it  would  be  natural  to 
develop  a  three  dimensional  version.  The  complexity  of  this  effort  is  not  trivial, 
even  with  the  two  dimensional  schemes  to  use  as  a  starting  points.  But  a  three 
dimensional  version  would  allow  computations  with  systems  that  do  not  possess 
azimuthal  symmetry  in  the  streaming  physics. 

Second,  some  small  modifications  to  FMP2DT  could  be  made.  This  could  include 
the  development  of  an  option  for  R0  geometry.  Modifications  to  the  input  scheme 
could  be  done  so  that  an  inhomogeneous  source  that’s  turned  on  and  off  would  not 
need  an  input  entry  describing  its  state  for  every  time  step.  At  present,  this  input 
is  very  lengthy.  Also,  it  would  be  desirable  to  include  an  option  in  FMP2DT  to 
C2.lcuiate  adjoint  fluxes.  These  and  several  other  minor  options  could  serve  as  some 
sort  of  academic  problem. 


106 


Appendix  A: 

Vector  and  Matrix  Definitions 


The  vector  Equations  (61)  and  (74)  are  the  result  of  applying  the  multigroup  ap¬ 
proximation  to  the  transport  equation.  They  have  been  defined  for  G  energy  groups. 
g  =  I  is  defined  as  the  highest  energy  group  and  g  =  G  is  defined  as  the  lowest, 
sometimes  referred  as  the  thermal  energy  group.  Therefore,  y  =  1, 2,  •  •  • ,  G.  We  now 
will  define  the  matrices  and  vectors  used  in  the  XZ  and  RZ  geometry  solutions.  All 
the  matrices  are  piece-wise  constant  in  an  interval. Ax  and  Az  or  Ar.  For  simplicity, 
their  following  definitions  will  not  carry  any  spatial  interval  subscripts. 

Define 


f  V>1  Vil-fl  1  1 

^t  ^.10  +  VI  At 

■^50 

_Vl-*2 

"so 

102  ’r'2-+2  1  1 

‘^50  +V2Af 

_y>2-*G 

"so 

.  vG  v^G-tG  1  1 

^so  +  voAt 

is  a  combination  of  Sf,  the  total  cross  section  for  the  energy  group, 
the  inscattering  cross  section  for  the  energy  group,  and  the  diagonal  entries 
come  from  combining  (addition)  to  the  V  matrix.  The  same  similar  combinations 
are  made  with  the  and  matrices. 

Define 


f  s;  -  2;,-'  + 


2”  =  •! 


At 
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2?  -  sr + 


V2  At 


-sfr' 


_2;-«  ...  +  J 


Define 


f  y>l  Vl-*1  A.  1 
^si  +  VI  At 

-sr 

^ 

-s:r’‘  s? 

y>2-»2  1  1 

+  „2  At 

-ssr® 

sf  j 

D'i  =  {Sii}~'and:D"  =  iDi' 

Define 


•  x'  S/ 

... 

X^v^Tfj  ' 

■■■ 

■■■ 

in  S'^  above  is- defined  as: 


X^  =  Xlil-^)  +  E)d0k  T-3 

k=l  X£. 


AfcAt, 


The  data  for  S°°,  S^°,  ^  and  S-^  are  given  in  a  cross  section  set  that  is  determined 

experimentally. 


Define  V: 


V=l 


vi  At 
0 


0 

1 


V2  At 


0 

0 


I  ^  .  ^t  ) 

Vg  is  the  neutron  velocity  related  to  the  neutron  energy.  At  is  the  time  discretization 
or  time  step  in  consideration. 
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For  the  vector  definitions,  their  spatial  superscript  i  j  will  refer  to  the  i  j*'^  mesh 


point.  The  superscript  (N  +  l)  will  relate  to  the  present  time  step  (At)  and  (iV)  will 
relate  to  the  previous  time  step.  The  energy  group  wdl  be  the  number  appearing 
under  the  time  step  designation. 


Define  the  ^qq  vectors: 


*00,y 


(W+1) 

V’ooly 

(W+1) 

^00;y 


'  (W) 

4o.y 

W 

*00.y  -  \ 


Define  the  ^ip  vectors: 


(AT+J) 


(W) 


^(AT+l)  _ 


(AT+I) 

i’lQij 

(W+l) 


*10,/  — 


(A'+l) 

^10% 


Define  the  vectors: 


(W) 


(N+l) 


r"'' 


(W) 

(N) 

Al/ 


(AT+l) 


Define  the  precursor  source: 


t,  l  +  Ai-At 


(W) 
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Then  the  precursor  vector  may  be  defined  as: 


Define  the  inhomogeneous  sources: 


(N) 

c'.-5 

(W) 


(N) 

pG 


(w+i) 
'JOO.y 
(N+l) 


(W+1) 
c  G 
■^OO.v,-  ) 


c(W+l)  _ 

°10,i  - 


(  (W+1) 

(N+1) 
O  2 


(A'+l) 

I  ^,0?,.  J 


r  (w+o  'v 


(JV+l) 

o  2 


(AT+l) 

Ui.?,-  J 


Note:  The  °  and  the  ^  matrices  were  presented  in  the  derivation  of  the  equations 
as  if  they  were  different.  Examination  of  the  entries  show  that  they  really  are  not. 
Thus,  the  diffusion  matrices,  and  are  also  the  same. 
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Appendix  B? 
Integration  Summaries 


The  linear  hat  functions  have  different  values  depending  on  the  interval  in  question. 
For  example,  for  an  interval  in  the  x  direction  from  z  —  1  to  z  +  1,  jB,(x)  has  two 
different  values.  This  can  be  seen  in  Figure  1.  Each  interior  spline  then  overlaps 
itself  and  its  two  nearest  neighbors.  Therefore,  thevB  splines  are  defined: 

In  the  X  direction  for  p  =  z  —  1,  z,  and:z  +  1: 


Rf_i(x)  = 

Xi  ■ 

-  X 

Xi- 

X;_l 

Bi{x)  = 

X  — 

X{-1 

X{- 

Xi-i 

Bi{x)  = 

®1+1 

—  X 

xi+1 

-Xi 

II 

+ 

X  - 

-Xi 

Sf+l 

—  Xi 

for  Xi-i  <Xi 
for  Xi^i  ^x  <Xi 
for  Xi  <x  <  x.’+i 
for  x;  <  a;  <  x,-+i 


In  the  z  direction  fov  q  =j  —  1,  j,  and  j  +  1: 


Z’  z 

Bj^-i{z)  =  — - -  for  Zj-i  <  z  <  Zj 

Zj  —  Zj-i 

Bj{z)  =  ~ — for  Zj^i  <  z  <  zj 
-  zj-i 

^  ^  ^  ^i+i 

^i+i  - 
Z  **  Z’ 

Bj^i{z)  = - —  for  Zj  <z  <  zj+i 

^i+i  - 


dR;_i(x) 

-1 

dx 

Xi  X,‘_i 

dBi{x) 

1 

dx 

Xi  X,'_i 

dBi{x) 

-1 

dx 

dJ5,+i(x) 

1 

dx 

dBj^i{z) 

-1 

dz 

Zj  -  Zj_i 

dBj{z)  _ 

1 

dz 

dBj{z) 

-1 

dz 

2j+l  -  Zj 

dBj^x{z) 

1 

dz 

""  ^i+i  -  Zj 

The  B  splines  in  the  z  direction  are  the  same  for  both  XZ  and  RZ  geometries. 


In  the  r  direction  for  p  =  ^  —  1,  i,  and  z  +  1 


B,(r)  = 

Bi(r)  = 


A«(r)  = 


r,-  —  r 


ri- 

r.-i 

r  — 

n-i 

r.-  - 

T’.-i 

—  r 

n-+i 

-  r,- 

r  - 

-r,- 

n+i 

-Ti 

for  r,-_i  <r  <Ti 


for  r,-_i  <r  <  ri 


for  r,-  <  r  <  r,+i 


for  r,-  <  r  <  r,+i 


_  -1 

dr  Ti  —  r,-_i 

__L_ 

dr  rj  —  r,-_i 

rf^.(r)  -1 

dr  r;+i  -  r; 

dBi^^jr)  _  1 

dr  ri+i  -  r; 


Now  define: 

=  re,-  —  rcf-i  and  hi  =  rcf+i  —  re,-  for  XZ  geometry. 

=  r,-  —  ri_i  and  /i,-  =  r,+i  —  r,-  for  RZ  geometry. 

—  Zj^i  and  hj  =  zj+i  —  Zj  for  both  geometries. 

For  integration  in  the  r  direction  with  an  r  factor  in  the  integrand: 
Define: 


Riii)  = 

R2(i)  = 

R3(i)  = 

RS)  = 
Rsd)  = 
Rc(i)  = 
Riii)  = 
Rs{i)  = 

Roii)  = 

^10  (i)  = 


i  (3r?  -  2r.r._i  -  rf_i) 
^  (r?+,  +  2rf+ir;  -  3r?) 
^  (^?+i  -  rf)  =Ri(i-  1) 


-1  I  Tj  +  r.-i 

2  \r,  -r,_i 

1  + 

2  Ir;  - r,_i / 


=  -Rsii) 


1  + 

2  Vr,+i  -  r{) 


-  (ri  +  2r,_i) 

-  {27-, •  +  7-;_i) 
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Bu(i)  =  (-^)  (ri+,+2r()  =  -i!9{i-l) 

=  (^)  (2n+,+r.)  =  -i?,„(i-l) 

Each  material  matrix  is  piece-wise  continuous  in  a  given  xz  or  rz  interval.  They 
possess  subscripts  identifying  the  specific  interval  where  they  belong.  Consider  the 
following  example.  Let  be  the  matrix  in  question.  Figure  32  shows  how  it  might 
fit  in  a  two  dimensional  slab  mesh  scheme.  Let  M^,  M^,  and  M'*  represent 


z 


Figure  32:  Piece- Wise  Material  Matrix  Example 
in  the  intervals  shown.  Then, 

M2  =  M**  =  S?P_i 

Thus,  is  continuous  in  the  interval  A,-_i  and  /ij_i,  is  continuous  in 

the  interval  /i;_i  and  hj,  is  continuous  in  the  interval  hi  and  hj,  and  SfjLj 
is  continuous  in  the  interval  /i,  and  hj-i.  All  other  material  matrices  arc  likewise 
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determined.  Now  the  integrations  are  addressed. 
In  the  X  direction: 


/  dBi^xix)  dBi{x)  , 

J  ~~di — 

1-1 

?  dBi{x)  dBi{x) 

J  ~di - dT''^ 

aTi-l 

ydBiix)dBiix) 

J  dx  dx 

Xi 

ydBi{x)dBi{x)^ 

J  —X - d^'^^ 


J  Bi-i  (x)  Bi{x)  dx 

xi-i 

xi 

j  Bi{x)  B{{x)  dx 

xi-t 

*i+j 

J  Bi{x)  Bt{x)  dx 

If 

if+i 

J  Bi+i{x)Bi{x)dx 


If 
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In  the  z  direction: 


J  dBi.,{z)  dBj^z) 

J  dz  dz 

r-i 

7  dBjjz)  dBj(z) 

J  dz  dz  " 

=3-1 

y  dBjjz)  dBjjz) 

J  dz  dz 

*3 

y  dBjjz)  dBjjz) 

J  dz  dz  ‘ 


•J 

‘/-I 

^3 

I  B,iz)Bi{z)dz 

=3-1 

=j*l 

J  ^siz)BMdz 
J  Bi:„(z)Bj{z)dz 

^3 


-J 
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In  the  r  direction: 


r; 

j  Bi(r)  r  dr 

fi-i 

r: 

j  Bi{r)  Bi{r)  r  dr 

••i-i 

n+i 

J  B{{r)  Bi{r)  r  dr 

r-i 

'■i+l 

j  Bi+i{r)Bi(r)rdr 

r; 


J 

ri-t 


r  dBj.ijr)  dBi{r) 


dr 


dr 


ri 

J 

ri-i 

J 

ri 

Tl+l 

I 

ri 

ri 


r,-i 


dBi{r)  dBi{r) 
dr  dr 

dB{{r)  dBijr) 
dr  dr 

dBi{r)  dBi{r) 
dr  dr 

dBi{r) 
dr 

dB{{r) 
dr 


J  Bi-.M 


/  £.(>•) 


r.--i 
r.+  i 


/  B,(r) 


r; 

’•t'+l 


/  Bwir) 


dBi(r) 

dr 

dr 


r  dr 


r  dr 


r  dr 


r  dr 


r  dr 


r  dr 


r  dr 


r  dr 


=  iJi(i) 
=  -BiCi) 

=  Mi) 
=  RS) 
=  Rsit) 

=  Mi) 

=  R7(i) 

=  Rsii) 
=  R,(i) 

=  i?io(0 
=  Rn{i) 
=  i?i2(i) 


ri 

Continued  on  the  next  page. 
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Some  of  the  integrands  in  the  radial  direction  did  not  have  a  r  factor  in  them 
because  it  was  canceled  out  by  a  ^  factor  in  the  same  integrand. 


1-1 

/  ^ 


iMlir  =  i 

dr  2 

dBiir) 


ri-t 

rj+i 


»■; 

*■1+1 


dr  =  ^ 
2 


dr  =  — — 
dr  2 


dr  =  — 
dr  2 


Finally  there  is  a  case  that  is  opposite  from  the  preceding  one. 


/ 

1-1 


dBi.tjr) 

dr 


Bi{r)  dr 


/ 


dBiir) 

dr 


Bi{r)  dr 


ri+i 

J 


r,- 


dBiir) 

dr 


Biir)  dr 


7‘  dBi+iir)  rjf 

J  — 


-1 

2 

1 

2 

-1 

2 

1 

2 
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Appendix  C: 

TWIGL  Cross  Section  Data 


The  following  is  a  summary  of  the  cross  sections  used  to  set  up  the  FEMP  calcula¬ 
tions  for  the  TWIGL  problem.  Consider  a  two  group  problem  with  a  Pi  approxi¬ 
mation.  The  pQ  transfer  matrix  is  defined  as: 

^00  _  wAt 

^2  V'2-»2  I  1 

~-^J0  iTAt  . 

And  the  Pi  transfer  matrix  is  defined  as: 


(3)  =  (3)  S/'  =  (3) 


s;  -  Eir' + 


vi  At 


-s?: 


51 


S?  - 


V2  J 


For  steady  state  calculations,  the  values  are  set  to  zero.  The  values  to  compute 
these  matrices  are  either  read  in  by  a  cross  section  tape,  tape  16,  or  read  in  the 
input  data  from  tape  5.  Consider  that  there  are  two  materials,  and  their  cross 
section  data  are  read  in  the  input  deck  (tape  5)  in  the  20**,  21**,  and  22**  arrays. 
Each  material  must  have  input  for  these  arrays.  For  the  TWIGL  problem,  the 
following  pages  show  the  input  for  the  steady  state  calculations.  For  the  FMP2DT 
calculations  at  time  /  =  0,  another  set  of  these  arrays  must  be  included  for  material 
3.  In  that  case,  material  3’s  arrays  would  be  exactly  like  material  I’s.  This  will 
change  for  t  >  0.005  sec  where  other  20**,  21**,  and  23**  arrays  must  be  entered, 
one  for  each  time  step,  reflecting  the  AS  for  each  entry  (as  described  in  the  text). 


Material  1: 


20**ARRAY: 

Sj  =  0.2398 
S?  =  0.6410 
I/S}  =  0.02555562 
I/S}  =  0.73259410 
S}  =  0.00 
S}  =  0.43 
X'  =  1-0 
x'  =  0.0 

21**ARRAY 

S}-*^  =  0.2148 
*0 

^  QQQgQ 
=  0.0000 
S}f2  =  0.2110 

22**ARRAY 

(3)S},-^  =  0.0 
(3)S}p2  =  0.0 
(3)S}p^  =  0.0 

(3)S}p2  =  0.0 
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Material  2: 


20**’ARRAY : 


= 

0.2801 

= 

0.4762 

I/S} 

= 

0.00511112 

z/S} 

= 

0.08518539 

= 

0.00 

= 

0.13 

= 

1.0 

0.0 

21**ARRAY 

=  0.2641 

50 

=  0.0060 
=  0.0000 
=  0.3462 

22**ARRAY 
(3)SJ7*'  =  0.0 

(3)SJr2  =  0.0 
(3)S2-i  =  0.0 
(3)S2f2  =  0.0 
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We  want  to  calculate  for  both  materials: 


= 


s;  -  Eir' 


-s?r'  : 
s?  -  . 


(Here  we  are  only  considering  steady  state  conditions.) 
Material  1: 


ol-^l 

^so 

=  0.2398- 

0.2148 

V>l-+2 

^40 

=  -  0.006 

V'2— ►! 

=  0.0 

V'2-»2 

^40 

=  0.6410- 

0.2110 

Thus: 


= 


0.0250  0.0000 
-0.006  0.4300 


Note,  it  is  no  accident  that  the  last  entry,  S?  —  =  0.4300  =  1 

matrix  for  a  Pi  approximation  is  defined  as: 


=  D^^=D={  (3) 


E?  - 


3(sj-sj;"^)  =  3(0.2398-0.0) 
3(S2-S2p2)  =  3(0.6410-0.0) 


Therefore; 


■  0.7194 

0.0000  ■ 

-1 

■  1.3900 

0.0000  ■ 

0.0000 

1.9230 

0.0000 

0.5200 

Iq.  The  diffusion 


-1 
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Material  2: 


Thus: 


SJ  -  = 


«0 


- 


-2' 


2-fl 

So 


0.2801  -  0.2641 
-  0.006 
0.0 

0.4762  -  0.3462 


S°°  = 


0.0160  0.0000 
-0.006  0.1300 


Note,  the  last  entry,  Sj  —  =  0.1300  =  S^.  The  diffusion  matrix  for  a  Pi 

approximation  is  defined  as: 


3(SJ-Sjp^)  =  3(0.2801  -0.0) 
3(S2-S2p2)  =  3(0.4762-0,0) 


Therefore; 


■  0.8403 

0.0000  ■ 

-1 

■  1.1901 

0.0000  ■ 

0.0000 

1.4286 

0.0000 

0.7000 
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The  MFP  values  yield  the  sensitivity  for  the  spatial  mesh  spacing  in  the  FEMP 
codes. 
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Material  1: 


Di  = 


1.0 


1.0 


3  3  (0.2398-0.000) 

Note,  this  is  the  D(l,l)  entry  in  the  diffusion  matrix. 


=  1.3900 


ii  = 


Di 


1.3900 


1.3900 


V  0.2398  -  0.2148  V  0.0250 
Note  the  removal  cross  section  for  group  1,  =  0.0250. 


=  7.4565 


M  —  — 


1.0 


D2  = 


SJ  0.2398 

1.0  1.0 


=  4.1701 


3  (S?  -  3  (0.6410  -  0.000) 

Note,  this  is  the  D(2,2)  entry  in  the  diffusion  matrix. 


=  0.5200 


i2  = 


D. 


0.5200 


0.5200 


S?  -  S2-2  y  0.6410  -  0.2110  V  0.4300 


=  1.0997 


Note  the  removal  cross  section  for  group  2,  2^  =  0.4300  =  2^. 


A2 


1.0  _  1.0 

2?  ~  0.6410 


1.5601 
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Material  2: 


I>i  = 


1.0 


1.0 


3  (S?  -  Si-i)  3  (0.2801  -  0.000) 
Wote,  this  is  the  D(l,l)  entry  in  the  diffusion  ma,trix. 


=  1.1901 


ii  = 


A 


1.1901 


1.1901 


^]  SJ  V  0-2801  -  0.2641  V  0.0160 

Note  the  removal  cross  section  for  group  1,  Sj  =  0.0160. 

1.0 


=  8.6243 


A  =H  = 

^  0.2801 

1.0 


=  3.5702 

1.0 


3  3  (0.4762  -  0.000) 

Note,  this  is  the  D  (2,2)  entry  in  the  diffusion  matrix. 


=  0.7000 


1 


Do 


0.7000 


S?  -  S2-2 


0.4762  -  0.3462 


fO.7000 

0.1300 


=  2.3205 


Note  the  removal  cross  section  for  group  2,  Sj  =  0.1300  = 


A2 


1.0  _  1.0 

S?  ”  0.4762 


2.1000 


Attached  are  FEMP2D  and  FMP2DT  input  decks  showing  the  20**,  21**,  and  22** 
array  input  for  the  TWIGL  problem. 
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FEMP2D  Sample  Input  Deck  (See  FEMP2D  Manual  for  details.) 


X-Z  Slab  Problem:  FEMP2D:  TWIGL  mesh  spacing. 

ISS 

TITLE 

1 

NGEOM 

3 

NOUTR 

0 

MADJ 

1 

LPN 

2 

NMAT 

2 

NNG 

0 

NPG 

1 

MPN 

0 

IHT 

0 

IHS 

0 

LTBL 

2 

MTL 

2 

MCRD 

0 

MANSN 

0 

MAMPX 

1 

NBYTE 

21 

NX 

4 

NY 

5 

NZONE 

0 

IB(1) 

0 

IB(2) 

0 

IB(3) 

0 

IB(4) 

2 

ISTRT 

4 

KSOLV 

500 

ITMXl 

200 

ITMX3 

0 

lACC 

0 

NPOW 

0 

NUPS 

0 

NS 

1 

IPX 

38 

NPOUT 

1 

IPFLX 

0 

2irk 

NRF 

l.OE-6 

EPS 

l.OE-5 

EPSK 

1.0 

XK 

8.37702E+16 

T 

SNORM 

loss 

1  2 

(Material  Numbers) 

llSS 

1  2 

(Nuclide  Numbers) 

12** 

1.0  1.0 

(Number  Dcn.sities) 

13SS 

1  2 

(Nuclide  IDs) 

14** 

(Prornpi  Chi  Spectrum) 

1.0 

X  X 

o;o 

17** 

(Group  Velocities) 

S.OE+06 

2.0E+05 

V2 

20** 

Material  1 

0.2398 

H 

0.6410 

2.555562E-02 

I/S> 

7.32S943E-01 

0.0 

i 

0.43 

1.0 
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0.0 

21** 

0.2148 

0.006 

0.0 

0.2110 

22** 

0.0 

0.0 

0.0 

0.0 

T 

20** 

0.2801 

0.4762 

5.11112E-03 

8.518539E-02 

0.0 

0.13 

1.0 

0.0 

21** 

0.2641 

0.006 

0.0 

0.3462 

22** 

0.0 

0.0 

0.0 

0.0 

T 

30** 

0.0  10.0  20.0  30.0  40.0  50.0 
60.0  70.0  80.0  90.0  100.0  110.0 
120.0  130.0  140.0  150.0  160.0  170.0 
180.0  190.0  200.0 
31** 

0.0  14.142  28.284  42.426 
3-lSS 

1112222333334444555 

1112222333334444555 

1112222333334444555 

35** 

FO.O 

36** 

10.0  20.0  30.0  40.0  50.0 
60.0  70.0  80.0  90.0  100.0  110.0 
120.0  130.0  140.0  150.0  160.0  170.0 
180.0  190.0 

10.0  20.0  30.0  40.0  50.0 

60.0  70.0  80.0  90.0  100.0  110.0 

120.0  130.0  140.0  150.0  160.0  170.0 

180.0  190.0 

37** 


14.142 

14.142 

14.142 

14.142 

14.142 

14.142 

14.142 

14.142 

14.142 

14.142 

14.142 

14.142 

14.142 

14.142 

14.142 

14.142 

14.142 

14.142 

14.142 

28.284 

28.284 

28.284 

28.284 

28.284 

28.284 

28.284 

28.284 

28.284 

28.281 

28.284 

28.284 

28.284 

28.284 

28.284 

28.281 

28.284 

28.284 

28.284 

T 


4 

SJo-*" 

mir^ 

(3)Sjp2 

(3)Eir‘ 

(3)Si-2 

Material  2 

i/Si 

i/Sj 

Si 

Si 

4 

^1—1 

SJo-* 

s5o-‘ 

s?r" 

(3)Ejp» 

(3)E1-2 

(3)s?r‘ 

(3)S5.-* 

X‘ Points 


Y  Points 


X  Print 


Y  Print 
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FMP2DT  Sample  Input  Deck  (See  FMP2DT  Manual  for  details.) 


X-Z  TWIGL  Problem:  (FMP2DT  run;  TVVIGL  mesh) 

TITLE 

ISS 

1 

NGEOM 

2 

NINT 

1 

ICOLD 

C 

LAM 

3 

NOUTR 

0 

MADJ 

1 

LPN 

3 

NMAT 

2 

NNG 

0 

NPG 

1 

MPN 

0 

IHT 

0 

IHS 

0 

LTBL 

3 

MTL 

3 

MCRD 

0 

MANSN 

0 

MAMPX 

1 

NBYTE 

21 

NX 

4 

NY 

5 

NZONE 

0 

IB(1) 

0 

IB(2) 

0 

IB(3) 

0 

IB(4) 

2 

ISTRT 

4 

KSOLV 

500 

ITMXl 

200 

ITMX3 

0 

lACC 

0 

NPOW 

C 

NUPS 

1 

IPX 

33 

NPOUT 

1 

IPFLX 

0 

NRF 

0 

IBAL 

2a* 

0.0 

TSTRT 

1.0E^ 

EPS 

1.0B-*I 

EPSK 

1.02C338 

XK 

8.37702E+1G 

SNORM 

T 

5$S 

(Material  Numbers) 

123 

6S$ 

(Nuclide  Numbers) 

123 

7** 

(Number  Densities) 

1.0  1.0  1.0 

sss 

(Nuclide  Ids) 

123 

0** 

(Prompt  Chi  Spectrum) 

1.0 

Xp 

0.0 

5 

xl 

12** 

(Group  Vdodties) 

5.0E+06 

ui 

2.0E+05 

n 

13** 

(Delayed  Neutron  Fractions) 

2.4700B-0'1 

ft 

1.38'15I>03 

ft 

1.2220&O3 

ft 
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2.6455I>03 

S3200E-04 

1.6900E^1 

0.0127 

0.0317 

0.1150 

0.3110 

1.4000 

3.8700 

15** 

1.0 


1.0 


1.0 


1.0 


1.0 

1.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

20** 

0.2398 

0.6410 

2.555562E-02 

7.325943&01 

0.0 

0.43 

1.0 

0.0 

21** 

0.2148 

0.006 

0.0 

0.2110 

22** 

0.0 

0.0 

0.0 

0.0 

T 

20** 

0.2^1 

0.4762 

5.11112E-03 

8.5I8539E-02 

0.0 

0.13 

1.0 

0.0 

21** 

0.2611 

0.006 

0.0 

0.3462 

22** 

0.0 

0.0 

0.0 

0.0 


ft 

ft 

ft 

(Decay  Constants) 

^2 

A* 

As 

As 

(Delayed  Chi  Spectnim) 

x\ 

xl 

4 

xl 

xi 

4 

4 

4 

4 

4 

4 

4 

1 

Vi 


Vl-I 
Vi-2 
V2— 1 
V2^2 


(3)sir* 

(3)sir’ 

(3)s;7-» 

(3)sir’ 


Material  2 
vl 


vi 


V2 


VI— I 


V2— I 
V^2 


(3)Eir> 


—2 
rJ— 1 


(3)Si. 


(3)Si- 

-  ‘T2i 

'•t 


C3)sjr" 
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T 

20** 

0.2398 

0.6410 

2.555562E-02 

7.325943E-01 

0.0 

0.43 

1.0 

0.0 

21** 

0.2148 

0.006 

0.0 

0.2110 

22** 

00 

0.0 

0.0 

0.0 

T 

30** 

0.0  10.0  20.0  30.0  40.0  50.0 
60.0  70.0  80.0  90.0  100.0  110.0 
120.0  130.0  140.0  150.0  160.0  170.0 
180.0  190.0  200.0 
31** 

0.0  14.142  28.284  42.426 
34$$ 

111  2-2  2233333444  4.5  55 
1  1  1  2  222333334444555 
1112222333334444555 
35** 

FO.O 

36** 

10.0  20.0  30.0  40.0  50.0 
60.0  70.0  80.0  90.0  100;0  110.0 
120.0  130.0  140.0  150.0  160.0  170.0 
180.0  190.0 

10.0  20.0  30.0  40.0  50.0 

60.0  70.0  80.0  90.0  100.0  110.0 

120.0  130.0  140.0  150.0  160.0  170.0 

180.0  190.0 

37** 


14.142 

14.142 

14.142 

14.142 

14.142 

14.142 

14.142 

14.142 

14.142 

14.142 

14.142 

14.142 

14.142 

14.142 

14.142 

14.142 

14.142 

14.142 

14.142 

28.284 

28.284 

28.284 

28.284 

28.284 

28.284 

28.284 

28.284 

28.284 

28.284 

28.284 

28.284 

28.284 

28.284 

28.284 

28.284 

28.284 

28.284 

28.284 

T 

40$$ 

0 

1 

0 

41** 

0.005 

T 

40$$ 

0 

5 

1 

41** 

0.010 


Material  3 

S| 

r/S? 

SJ 

SI 

Xp 

SJo-*‘ 

SV‘ 

s?„-^ 

(3)sjr‘ 

(3)sjr* 

(3)s?r‘ 

(3)s?r^ 

X  Points 


Y  Points 


X  Print 


Y  Print 


INT=1 

NS 

NTIM 

IXSEC 

TFIN 


INT=2 

NS 

NTIM 

IXSEC 

TFIN 
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T 

20Mc 

0.0 

0.0 

0.0 

0.0 

0.0 

-0.002 

0.0 

0.0 

21** 

0.0 

0.0 

0.0 

0.002 

22** 

FO.O 

T 

20** 

FO.O 

21** 

FO.O 

22** 

FO.O 

T 

20** 

FO.O 

21** 

FO.O 

22** 

FO.O 

T 

20** 

0.0 

0.0 

0.0 

0.0 

0.0 

-0.002 

0.0 

0.0 

21** 

0.0 

0.0 

0.0 

0.002 

22** 

FO.O 

T 

20** 

FO.O 

21** 

FO.O 

22** 

FO.O 

T 

20** 

FO.O 

21** 

FO.O 

22** 

FO.O 

T 

20** 

0.0 

0.0 


Material  1  (ITT=l) 

§ 

Si 

Si 

Xp 


sjr* 

SJo-" 

S?o-‘ 

S?o“^ 


Material  2  (ITT=1) 


Material  3  (ITT=l) 


Material  1  (ITT=2) 
j/S‘, 


«/S' 

Si' 

Si 

Xp 


i 


So- 

Sio-* 

si„-** 


Material  2  (ITT=2) 


Material  3  (ITT=2) 


Material  1  (1TT=3) 
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0.0 

vyA 

0.0 

0.0 

Si 

-0.002 

Si 

0.0 

Xp 

0.0 

21** 

0.0 

0.0 

SJo-* 

0.0 

sio"*' 

0.002 

s?„-^ 

22** 

FO.O 

T 

20** 

Material  2  (ITT=3) 

FO.O 

21** 

FO.O 

22** 

FO.O 

T 

20** 

Material  3  (ITT=3) 

FO.O 

21** 

FO.O 

22** 

FO.O 

f' 

20** 

Material  1  (rTT=4) 

0.0 

§ 

0.0 

0.0 

i/S* 

pS? 

0.0 

0.0 

-0.002 

0.0 

0.0 

21** 

i 

i 

0.0 

s5o“‘ 

0,0 

sJo-*" 

0.0 

s?o-‘ 

0.002 

22** 

FO.O 

T 

20** 

Material  2  (ITT=4) 

FO.O 

21** 

FO.O 

22** 

FO.O 

T 

20** 

Material  3  (ITT=4) 

FO.O 

21** 

FO.O 

22** 

FO.O 

T 

20** 

Material  1  (ITT=5) 

0.0 

i/Sl 

u-A 

0.0 

0.0 

0.0 

0.0 

Si' 

-0.002 

Si 
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0.0 

Xp 

0.0 

21:^ 

0.0 

£1-1 

0.0 

Sio-* 

0.6- 

0.002 

22** 

FO.O 

T 

20** 

Material  2  (ITT=5) 

FO.O 

21** 

FO.O 

22** 

FO.O 

T 

20** 

Material  3  (ITT=5) 

FO.O 

21** 

FO.O 

22** 

FO.O 

T 

Note:  ITT=TIME  INTERVAL;  ITT=1,*  •  •,  NINT 

INT=TIME  STEP  IN  THAT  INTERVAL;  INT=1,.  •  •,  NTIM 
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